This notebook looks at the effect of changing different parameters on the tides.
Runs:
Measured amplitude/phase from Foreman's Discovery Islands and 2004 paper are included.
Complex differences are from the Foreman inversion method in 2004 paper.
Uses the same curve fitting technique that Susan wrote.
# 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
import collections
# grid
bathy, X, Y = tidetools.get_SS_bathy_data()
currently comparing tide_slip with the control case tide_flux_M2K1.
name1 = '/data/nsoontie/MEOPAR/SalishSea/results/tides/CBase/'
name2 = '/data/nsoontie/MEOPAR/SalishSea/results/tides/tide_flux_M2K1/'
# M2
M2freq = 28.984106 # degrees per hour
M2freq = M2freq*np.pi/180. # radians per hour
K1freq = 15.041069*np.pi/180.
location = ("PortRenfrew","Victoria","FulfordHarbour","PortTownsend", "PointAtkinson",
"PowellRiver", "Lund","CampbellRiver","YorkeIsland")
numsta=9
# initialize fit
M2amp = 1.; M2pha = 180.
K1amp = 1.; K1pha = 180.
# 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.))
# 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
#measurements from Foreman
#note some inconsistencies with Campbell River. Different measurements in different papers
#saved in /data/nsoontie/MEOPAR/analysis/compare_tides/Foreman_obs.csv
M2_amp = {"PortRenfrew": .708, "Victoria": .37, "PointAtkinson": .919, "PortTownsend": .684,
"FulfordHarbour": .582,
"CampbellRiver":.825, "YorkeIsland": 1.171, "PowellRiver": 1.007, "Lund": 1.022}
M2_pha = {"PortRenfrew": 241.1, "Victoria": 317.7, "PointAtkinson": 31, "PortTownsend": 350.5,
"FulfordHarbour": 12.7,
"CampbellRiver":18.4, "YorkeIsland": 271.8, "PowellRiver": 34.3, "Lund": 35.4}
K1_amp = {"PortRenfrew": .453, "Victoria": .636, "PointAtkinson": .862, "PortTownsend": .764,
"FulfordHarbour": .753,
"CampbellRiver":.846, "YorkeIsland": .558, "PowellRiver": .904, "Lund": .889}
K1_pha = {"PortRenfrew": 254.1, "Victoria": 269.6, "PointAtkinson": 286.2, "PortTownsend": 270.8,
"FulfordHarbour": 280,
"CampbellRiver":284, "YorkeIsland": 260, "PowellRiver": 286.6, "Lund":287.9}
ts = 240 # last half of serie
data={}
data[name1]={'M2 Amp':{}, 'M2 Pha':{}, 'K1 Amp':{}, 'K1 Pha':{}}
data[name2]={'M2 Amp':{}, 'M2 Pha':{}, 'K1 Amp':{}, 'K1 Pha':{}}
for name in (name1,name2):
plt.figure(figsize=(12,9))
print name
for stn in range(numsta):
fT1 = NC.Dataset(name+location[stn]+'.nc','r')
time = fT1.variables["time_counter"][:]/3600. # want hours not seconds
ssh = fT1.variables["sossheig"][:,0,0]
fitted, cov = curve_fit(double,time[ts:],ssh[ts:])
if fitted[0] < 0:
fitted[0] = -fitted[0]
fitted[1] = fitted[1]+180.
print location[stn]
print " M2 original {0:.3f} m and {1:.0f} degrees.".format(fitted[0], fitted[1])
print " M2 corrected {0:.3f} m and {1:.0f} degrees.".format(fitted[0]*M2ft, fitted[1]+M2uvt)
print " M2 measured {0:.3f} m and {1:.0f} degrees.".format(M2_amp[location[stn]], M2_pha[location[stn]])
data[name]['M2 Amp'][location[stn]] = fitted[0]*M2ft
data[name]['M2 Pha'][location[stn]] = fitted[1]+M2uvt
if data[name]['M2 Pha'][location[stn]] > 360:
data[name]['M2 Pha'][location[stn]] = data[name]['M2 Pha'][location[stn]]-360
print " K1 original {0:.3f} m and {1:.0f} degrees.".format(fitted[2],fitted[3])
print " K1 corrected {0:.3f} m and {1:.0f} degrees.".format(fitted[2]*K1ft, fitted[3]+K1uvt)
print " K1 measured {0:.3f} m and {1:.0f} degrees.".format(K1_amp[location[stn]], K1_pha[location[stn]])
data[name]['K1 Amp'][location[stn]] = fitted[2]*K1ft
data[name]['K1 Pha'][location[stn]] = fitted[3]+K1uvt
if data[name]['K1 Pha'][location[stn]] > 360:
data[name]['K1 Pha'][location[stn]] = data[name]['K1 Pha'][location[stn]]-360
plt.subplot(3,3,stn+1)
plt.plot (time, ssh, time, double(time,fitted[0], fitted[1], fitted[2], fitted[3]))
plt.title(location[stn])
/data/nsoontie/MEOPAR/SalishSea/results/tides/CBase/ PortRenfrew M2 original 0.713 m and -10 degrees. M2 corrected 0.704 m and 236 degrees. M2 measured 0.708 m and 241 degrees. K1 original 0.406 m and -43 degrees. K1 corrected 0.426 m and 254 degrees. K1 measured 0.453 m and 254 degrees. Victoria M2 original 0.271 m and 66 degrees. M2 corrected 0.268 m and 312 degrees. M2 measured 0.370 m and 318 degrees. K1 original 0.575 m and -25 degrees. K1 corrected 0.604 m and 272 degrees. K1 measured 0.636 m and 270 degrees. FulfordHarbour M2 original 0.463 m and 133 degrees. M2 corrected 0.457 m and 379 degrees. M2 measured 0.582 m and 13 degrees. K1 original 0.671 m and -14 degrees. K1 corrected 0.704 m and 283 degrees. K1 measured 0.753 m and 280 degrees. PortTownsend M2 original 0.544 m and 100 degrees. M2 corrected 0.537 m and 346 degrees. M2 measured 0.684 m and 350 degrees. K1 original 0.678 m and -25 degrees. K1 corrected 0.712 m and 272 degrees. K1 measured 0.764 m and 271 degrees. PointAtkinson M2 original 0.876 m and 150 degrees. M2 corrected 0.866 m and 396 degrees. M2 measured 0.919 m and 31 degrees. K1 original 0.783 m and -7 degrees. K1 corrected 0.822 m and 289 degrees. K1 measured 0.862 m and 286 degrees. PowellRiver M2 original 0.948 m and 152 degrees. M2 corrected 0.937 m and 398 degrees. M2 measured 1.007 m and 34 degrees. K1 original 0.806 m and -6 degrees. K1 corrected 0.847 m and 290 degrees. K1 measured 0.904 m and 287 degrees. Lund M2 original 0.972 m and 153 degrees. M2 corrected 0.961 m and 398 degrees. M2 measured 1.022 m and 35 degrees. K1 original 0.812 m and -6 degrees. K1 corrected 0.854 m and 290 degrees. K1 measured 0.889 m and 288 degrees. CampbellRiver M2 original 0.719 m and 130 degrees. M2 corrected 0.710 m and 376 degrees. M2 measured 0.825 m and 18 degrees. K1 original 0.772 m and -12 degrees. K1 corrected 0.811 m and 284 degrees. K1 measured 0.846 m and 284 degrees. YorkeIsland M2 original 1.651 m and 25 degrees. M2 corrected 1.631 m and 271 degrees. M2 measured 1.171 m and 272 degrees. K1 original 0.628 m and -52 degrees. K1 corrected 0.660 m and 244 degrees. K1 measured 0.558 m and 260 degrees. /data/nsoontie/MEOPAR/SalishSea/results/tides/tide_flux_M2K1/ PortRenfrew M2 original 0.724 m and -3 degrees. M2 corrected 0.716 m and 243 degrees. M2 measured 0.708 m and 241 degrees. K1 original 0.461 m and -39 degrees. K1 corrected 0.484 m and 258 degrees. K1 measured 0.453 m and 254 degrees. Victoria M2 original 0.294 m and 67 degrees. M2 corrected 0.291 m and 313 degrees. M2 measured 0.370 m and 318 degrees. K1 original 0.651 m and -20 degrees. K1 corrected 0.684 m and 276 degrees. K1 measured 0.636 m and 270 degrees. FulfordHarbour M2 original 0.445 m and 135 degrees. M2 corrected 0.440 m and 381 degrees. M2 measured 0.582 m and 13 degrees. K1 original 0.761 m and -9 degrees. K1 corrected 0.799 m and 288 degrees. K1 measured 0.753 m and 280 degrees. PortTownsend M2 original 0.548 m and 103 degrees. M2 corrected 0.542 m and 348 degrees. M2 measured 0.684 m and 350 degrees. K1 original 0.770 m and -20 degrees. K1 corrected 0.809 m and 276 degrees. K1 measured 0.764 m and 271 degrees. PointAtkinson M2 original 0.847 m and 156 degrees. M2 corrected 0.837 m and 401 degrees. M2 measured 0.919 m and 31 degrees. K1 original 0.892 m and -2 degrees. K1 corrected 0.937 m and 295 degrees. K1 measured 0.862 m and 286 degrees. PowellRiver M2 original 0.921 m and 157 degrees. M2 corrected 0.910 m and 403 degrees. M2 measured 1.007 m and 34 degrees. K1 original 0.918 m and -1 degrees. K1 corrected 0.965 m and 295 degrees. K1 measured 0.904 m and 287 degrees. Lund M2 original 0.947 m and 158 degrees. M2 corrected 0.935 m and 404 degrees. M2 measured 1.022 m and 35 degrees. K1 original 0.926 m and -1 degrees. K1 corrected 0.972 m and 295 degrees. K1 measured 0.889 m and 288 degrees. CampbellRiver M2 original 0.668 m and 142 degrees. M2 corrected 0.660 m and 388 degrees. M2 measured 0.825 m and 18 degrees. K1 original 0.864 m and -5 degrees. K1 corrected 0.908 m and 291 degrees. K1 measured 0.846 m and 284 degrees. YorkeIsland M2 original 1.315 m and 22 degrees. M2 corrected 1.299 m and 268 degrees. M2 measured 1.171 m and 272 degrees. K1 original 0.605 m and -42 degrees. K1 corrected 0.635 m and 254 degrees. K1 measured 0.558 m and 260 degrees.
plt.figure(figsize=(8,10))
plt.pcolormesh(X,Y,bathy)
for stn in range(numsta):
fT1 = NC.Dataset(name1+ '/' +location[stn]+'.nc','r')
lat = fT1.variables["nav_lat"]
lon = fT1.variables["nav_lon"]
plt.plot(lon[0,0],lat[0,0],'o',label=location[stn])
plt.legend(bbox_to_anchor=(1.05, 1), loc=2)
<matplotlib.legend.Legend at 0x7792710>
print 'Comparison of model with observations. New run and old run'
print 'new case / control case'
for stn in range(numsta):
print location[stn]
print 'M2 Amplitude ratios', data[name1]['M2 Amp'] [location[stn]]/M2_amp[location[stn]],'/', data[name2]['M2 Amp'] [location[stn]]/M2_amp[location[stn]]
print 'M2 Phase difference', data[name1]['M2 Pha'] [location[stn]]-M2_pha[location[stn]],'/', data[name2]['M2 Pha'] [location[stn]]-M2_pha[location[stn]]
print 'K1 Amplitude ratios', data[name1]['K1 Amp'] [location[stn]]/K1_amp[location[stn]],'/', data[name2]['K1 Amp'] [location[stn]]/K1_amp[location[stn]]
print 'K1 Phase difference', data[name1]['K1 Pha'] [location[stn]]-K1_pha[location[stn]],'/', data[name2]['K1 Pha'] [location[stn]]-K1_pha[location[stn]]
Comparison of model with observations. New run and old run new case / control case PortRenfrew M2 Amplitude ratios 0.994675505771 / 1.01067243011 M2 Phase difference -5.22977035859 / 1.73727327462 K1 Amplitude ratios 0.94052992584 / 1.06844072961 K1 Phase difference -0.397732087524 / 3.5095702148 Victoria M2 Amplitude ratios 0.724381637118 / 0.785786079038 M2 Phase difference -6.18040556011 / -4.58577956585 K1 Amplitude ratios 0.949737443569 / 1.0753123651 K1 Phase difference 2.13513025896 / 6.70516122422 FulfordHarbour M2 Amplitude ratios 0.78560120927 / 0.755893196611 M2 Phase difference 5.87398127897 / 8.68375165347 K1 Amplitude ratios 0.935500987594 / 1.06133945812 K1 Phase difference 2.65744894027 / 7.72640833731 PortTownsend M2 Amplitude ratios 0.785588007637 / 0.791798921169 M2 Phase difference -4.76963831667 / -2.0920151292 K1 Amplitude ratios 0.932479302751 / 1.05914281174 K1 Phase difference 0.783079183412 / 5.64023580378 PointAtkinson M2 Amplitude ratios 0.941983281501 / 0.910723977477 M2 Phase difference 5.03941319343 / 10.4498037099 K1 Amplitude ratios 0.954133140539 / 1.08733910595 K1 Phase difference 3.07865388032 / 8.45871089004 PowellRiver M2 Amplitude ratios 0.930265314182 / 0.903389185881 M2 Phase difference 3.39847318335 / 9.01985356868 K1 Amplitude ratios 0.936636579268 / 1.06722105368 K1 Phase difference 3.23051699595 / 8.69156746615 Lund M2 Amplitude ratios 0.939879647077 / 0.915009955973 M2 Phase difference 3.05836519504 / 8.77922116971 K1 Amplitude ratios 0.960096127659 / 1.09380330299 K1 Phase difference 2.05317903978 / 7.54613816065 CampbellRiver M2 Amplitude ratios 0.86075257339 / 0.800293338556 M2 Phase difference -2.63984424578 / 9.65507386849 K1 Amplitude ratios 0.958922121347 / 1.07299259656 K1 Phase difference 0.215662433559 / 7.29773001041 YorkeIsland M2 Amplitude ratios 1.39295282217 / 1.10946912206 M2 Phase difference -0.622717513696 / -3.91123120463 K1 Amplitude ratios 1.18288064357 / 1.13830066451 K1 Phase difference -15.6641741972 / -5.77513949076
Above: If phase change is above 360 then just subtract 360.
plt.figure(figsize=(10,8))
#There is a more efficient way to do this. My Dict objects aren't keeping the order.
obs=[]; one=[]; two=[]
for stn in range(numsta):
obs.append(M2_amp[location[stn]])
one.append(data[name1]['M2 Amp'][location[stn]])
two.append(data[name2]['M2 Amp'][location[stn]])
plt.subplot(2,2,1)
plt.plot(one, '-bo', label = 'new')
plt.plot(two, 'g-o', label = 'old')
plt.plot(obs, 'r-o', label = 'measured')
plt.title('M2 Amplitude')
plt.subplot(2,2,2)
obs=[]; one=[]; two=[]
for stn in range(numsta):
obs.append(K1_amp[location[stn]])
one.append(data[name1]['K1 Amp'][location[stn]])
two.append(data[name2]['K1 Amp'][location[stn]])
plt.plot(one, '-bo', label = 'new')
plt.plot(two, 'g-o', label = 'old')
plt.plot(obs, 'r-o', label = 'measured')
plt.title('K1 Amplitude')
plt.subplot(2,2,3)
obs=[]; one=[]; two=[]
for stn in range(numsta):
obs.append(M2_pha[location[stn]])
one.append(data[name1]['M2 Pha'][location[stn]])
two.append(data[name2]['M2 Pha'][location[stn]])
plt.plot(one, '-bo', label = 'new')
plt.plot(two, 'g-o', label = 'old')
plt.plot(obs, 'r-o', label = 'measured')
plt.title('M2 Phase')
plt.subplot(2,2,4)
obs=[]; one=[]; two=[]
for stn in range(numsta):
obs.append(K1_pha[location[stn]])
one.append(data[name1]['K1 Pha'][location[stn]])
two.append(data[name2]['K1 Pha'][location[stn]])
plt.plot(one, '-bo', label = 'new')
plt.plot(two, 'g-o', label = 'old')
plt.plot(obs, 'r-o', label = 'measured')
plt.title('K1 Phase')
for stn in range(numsta):
print stn, location[stn]
0 PortRenfrew 1 Victoria 2 FulfordHarbour 3 PortTownsend 4 PointAtkinson 5 PowellRiver 6 Lund 7 CampbellRiver 8 YorkeIsland
blue = newest run, green = older, red= measurements
** Summary **
Flux increase at west (not shown in this analysis)
** Bottom Friction** - 3e-3
** Viscosity Lower \nu=15 **
** Bottom Friction** - 1e-3
** K1_phase2 **
** K1 Amp**
** M2 Phase**
** tide_slip **
** tide_slipH**
** tide_bf0 **
** K1 amp/ phase at Port Renfrew **
Complex differences: Comparison for this run
Note: Think about a way to make it easier to compare between all runs?
def complex_diff(Ao,go,Am,gm):
#calculates complex differences between observations and model
import numpy as np
import math
D = np.sqrt((Ao*np.cos(math.radians(go))-Am*np.cos(math.radians(gm)))**2 +
(Ao*np.sin(math.radians(go))-Am*np.sin(math.radians(gm)))**2)
return D
Foreman's complex differences.
M2: From Foreman et al 2004 inversion method (Table 3, column 2)
K1: From Foreman et al 2004 (Table 1). The inversion analysis was only applied to M2.
Northern tides (Lund, Powell River is from Foreman et al Discovery Island, 2012)
Foreman_M2= {"PortRenfrew": .007, "Victoria": .003, "PointAtkinson": .005, "PortTownsend": .009,
"FulfordHarbour": .011,
"CampbellRiver":.142, "YorkeIsland": .03, "PowellRiver": .01, "Lund": .014}
Foreman_K1= {"PortRenfrew": .031, "Victoria": .053, "PointAtkinson": .207, "PortTownsend": .062,
"FulfordHarbour": .128,
"CampbellRiver":.197, "YorkeIsland": .097, "PowellRiver": .033, "Lund": .008}
plt.figure(figsize=(10,6))
D_M2=[]
D_K1=[]
name=name1
for stn in range(numsta):
D_M2.append(complex_diff(M2_amp[location[stn]],M2_pha[location[stn]]
,data[name]['M2 Amp'][location[stn]], data[name]['M2 Pha'][location[stn]] ) )
D_K1.append(complex_diff(K1_amp[location[stn]],K1_pha[location[stn]]
,data[name]['K1 Amp'][location[stn]], data[name]['K1 Pha'][location[stn]] ) )
plt.subplot(1,2,1)
plt.plot(D_M2,'ob-')
plt.title('M2 Complex Differences')
plt.subplot(1,2,2)
plt.plot(D_K1,'ob-',label='new run')
plt.title('K1 Complex Differences')
D_M2=[]
D_K1=[]
F_M2=[]; F_K1=[]
name=name2
for stn in range(numsta):
D_M2.append(complex_diff(M2_amp[location[stn]],M2_pha[location[stn]]
,data[name]['M2 Amp'][location[stn]], data[name]['M2 Pha'][location[stn]] ))
D_K1.append(complex_diff(K1_amp[location[stn]],K1_pha[location[stn]]
,data[name]['K1 Amp'][location[stn]], data[name]['K1 Pha'][location[stn]] ))
F_M2.append(Foreman_M2[location[stn]])
F_K1.append(Foreman_K1[location[stn]])
plt.subplot(1,2,1)
plt.plot(D_M2,'og-')
plt.plot(F_M2,'or-')
plt.title('M2 Complex Differences')
plt.subplot(1,2,2)
plt.plot(D_K1,'og-',label='old run')
plt.plot(F_K1,'or-',label='Foreman')
plt.title('K1 Complex Differences')
plt.legend(loc=2)
print 'blue:', name1
print 'green:', name2
print 'red:', 'Foreman'
print 'Stations:'
for stn in range(numsta):
print stn, location[stn]
blue: /data/nsoontie/MEOPAR/SalishSea/results/tides/CBase/ green: /data/nsoontie/MEOPAR/SalishSea/results/tides/tide_flux_M2K1/ red: Foreman Stations: 0 PortRenfrew 1 Victoria 2 FulfordHarbour 3 PortTownsend 4 PointAtkinson 5 PowellRiver 6 Lund 7 CampbellRiver 8 YorkeIsland
Note: Foreman's Discovery Island model gives much better results for Yorke Island, Campbell River...
Foreman has better commplex differences in general. This is the inversion technique. I could compare with Foreman et al 1995 but there aren't as many stations.
K1 is variable. Most of the time Foreman's are better. Remember these are from the tide3d code in Foreman et al 2004 (without inversion correction).