Study Tides: Multi-constituents : Fit to Tide Frequencies
In [2]:
# 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
In [3]:
# grid
bathy, X, Y = tidetools.get_SS_bathy_data()
In [27]:
name1 = '/ocean/sallen/allen/research/Meopar/myResults/SalishSea/M2K1only/'
name2 = '/ocean/nsoontie/MEOPAR/SalishSea/results/tide_flux_M2K1/'
# M2
M2freq = 28.984106 # degrees per hour
M2freq = M2freq*np.pi/180. # radians per hour
K1freq = 15.041069*np.pi/180.
print M2freq, K1freq
location = ("PortRenfrew","Victoria","Sandheads","PowellRiver","YorkeIsland")
0.505868080447 0.26251617707
Now some useful functions: first fit, then phase
In [5]:
# initialize fit
M2amp = 1.; M2pha = 180.
K1amp = 1.; K1pha = 180.
m = 0.; b = 0
In [6]:
# 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 [7]:
# function for fit, with linear part too
def doubleplus(x, M2amp, M2pha, K1amp, K1pha, m, b):
return (M2amp*np.cos(M2freq*x-M2pha*np.pi/180.)+
K1amp*np.cos(K1freq*x-K1pha*np.pi/180.) + m*x + b)
In [28]:
# function for fit, with M4
def doublenl(x, M2amp, M2pha, K1amp, K1pha, M4amp, M4pha):
return (M2amp*np.cos(M2freq*x-M2pha*np.pi/180.)+
K1amp*np.cos(K1freq*x-K1pha*np.pi/180.)+
M4amp*np.cos(2*M2freq*x-M4pha*np.pi/180.))
In [10]:
# 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
In [12]:
plt.figure(figsize=(12,9))
for stn in range(5):
fT1 = NC.Dataset(name1+location[stn]+'.nc','r')
time1 = fT1.variables["time_counter"][:]/3600. # want hours not seconds
ssh1 = fT1.variables["sossheig"][:,0,0]
fT2 = NC.Dataset(name2+location[stn]+'.nc','r')
time2 = fT2.variables["time_counter"][:]/3600. # want hours not seconds
ssh2 = fT2.variables["sossheig"][:,0,0]
plt.subplot(2,3,stn+1)
plt.plot (time1, ssh1, time2, ssh2)
plt.title(location[stn])
Comparing the results with and without the changes in the north: only difference are in the north and major difference is decrease in M2 tide magnitude
In [15]:
ts = 240 # last half of series
name = name1
plt.figure(figsize=(12,9))
for stn in range(5):
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 " 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)
plt.subplot(2,3,stn+1)
plt.plot (time, ssh, time, double(time,fitted[0], fitted[1], fitted[2], fitted[3]))
plt.title(location[stn])
PortRenfrew
M2 original 0.718 m and -3 degrees.
M2 corrected 0.709 m and 243 degrees.
K1 original 0.462 m and -39 degrees.
K1 corrected 0.485 m and 258 degrees.
Victoria
M2 original 0.292 m and 69 degrees.
M2 corrected 0.289 m and 315 degrees.
K1 original 0.652 m and -20 degrees.
K1 corrected 0.685 m and 276 degrees.
Sandheads
M2 original 0.821 m and 154 degrees.
M2 corrected 0.811 m and 400 degrees.
K1 original 0.883 m and -2 degrees.
K1 corrected 0.928 m and 295 degrees.
PowellRiver
M2 original 0.927 m and 157 degrees.
M2 corrected 0.916 m and 403 degrees.
K1 original 0.919 m and -1 degrees.
K1 corrected 0.965 m and 295 degrees.
YorkeIsland
M2 original 1.667 m and 25 degrees.
M2 corrected 1.647 m and 271 degrees.
K1 original 0.627 m and -50 degrees.
K1 corrected 0.659 m and 246 degrees.
Check if we should be adding a linear term and a constant to fit¶
In [16]:
ts = 240 # last half of series
name = name1
plt.figure(figsize=(12,9))
for stn in range(5):
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.
fitted2, cov2 = curve_fit(doubleplus,time[ts:],ssh[ts:])
if fitted2[0] < 0:
fitted2[0] = -fitted2[0]
fitted2[1] = fitted2[1]+180.
print location[stn]
print " M2 corrected {0:.3f} m and {1:.0f} degrees.".format(fitted[0]*M2ft, fitted[1]+M2uvt)
print " M2 better? {0:.3f} m and {1:.0f} degrees.".format(fitted2[0]*M2ft, fitted2[1]+M2uvt)
print " K1 corrected {0:.3f} m and {1:.0f} degrees.".format(fitted[2]*K1ft, fitted[3]+K1uvt)
print " K1 better? {0:.3f} m and {1:.0f} degrees.".format(fitted2[2]*K1ft, fitted2[3]+K1uvt)
plt.subplot(2,3,stn+1)
plt.plot (time, ssh, time, double(time,fitted[0], fitted[1], fitted[2], fitted[3]),
time, doubleplus(time,fitted2[0],fitted2[1],fitted2[2],fitted2[3],fitted2[4],fitted2[5]))
plt.title(location[stn])
PortRenfrew
M2 corrected 0.709 m and 243 degrees.
M2 better? 0.709 m and 243 degrees.
K1 corrected 0.485 m and 258 degrees.
K1 better? 0.481 m and 258 degrees.
Victoria
M2 corrected 0.289 m and 315 degrees.
M2 better? 0.290 m and 315 degrees.
K1 corrected 0.685 m and 276 degrees.
K1 better? 0.682 m and 277 degrees.
Sandheads
M2 corrected 0.811 m and 400 degrees.
M2 better? 0.811 m and 400 degrees.
K1 corrected 0.928 m and 295 degrees.
K1 better? 0.927 m and 294 degrees.
PowellRiver
M2 corrected 0.916 m and 403 degrees.
M2 better? 0.917 m and 403 degrees.
K1 corrected 0.965 m and 295 degrees.
K1 better? 0.965 m and 295 degrees.
YorkeIsland
M2 corrected 1.647 m and 271 degrees.
M2 better? 1.647 m and 271 degrees.
K1 corrected 0.659 m and 246 degrees.
K1 better? 0.659 m and 246 degrees.
Adding a linear and constant term make a maximum of 4 mm and/or 1 degree difference
What difference does length of time series make?¶
In [17]:
for ts in range(96,288,24):
name = name1
stn = 2
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 " M2 corrected {0:.3f} m and {1:.0f} degrees.".format(fitted[0]*M2ft, fitted[1]+M2uvt), " K1 corrected {0:.3f} m and {1:.0f} degrees.".format(fitted[2]*K1ft, fitted[3]+K1uvt);
M2 corrected 0.812 m and 400 degrees. K1 corrected 0.930 m and 294 degrees.
M2 corrected 0.813 m and 400 degrees. K1 corrected 0.935 m and 294 degrees.
M2 corrected 0.812 m and 400 degrees. K1 corrected 0.930 m and 295 degrees.
M2 corrected 0.812 m and 400 degrees. K1 corrected 0.924 m and 294 degrees.
M2 corrected 0.811 m and 400 degrees. K1 corrected 0.928 m and 294 degrees.
M2 corrected 0.814 m and 400 degrees. K1 corrected 0.934 m and 294 degrees.
M2 corrected 0.811 m and 400 degrees. K1 corrected 0.928 m and 295 degrees.
M2 corrected 0.813 m and 400 degrees. K1 corrected 0.923 m and 294 degrees.
As we go from analyzing 4 days to 2 days (last 4 to last 2) the variation in the amplitude of M2 is maximum 3 mm and no change in phase. The variation in the amplitude of K1 is 11 mm (1 cm) and 1 degree in phase.
Do we need M4 at, say, Victoria?¶
In [31]:
ts = 240 # last half of series
name = name1
plt.figure(figsize=(5,5))
stn = 1
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.
fitted2, cov2 = curve_fit(doublenl,time[ts:],ssh[ts:])
if fitted2[0] < 0:
fitted2[0] = -fitted2[0]
fitted2[1] = fitted2[1]+180.
print location[stn]
print " M2 corrected {0:.3f} m and {1:.0f} degrees.".format(fitted[0]*M2ft, fitted[1]+M2uvt)
print " M2 better? {0:.3f} m and {1:.0f} degrees.".format(fitted2[0]*M2ft, fitted2[1]+M2uvt)
print " K1 corrected {0:.3f} m and {1:.0f} degrees.".format(fitted[2]*K1ft, fitted[3]+K1uvt)
print " K1 better? {0:.3f} m and {1:.0f} degrees.".format(fitted2[2]*K1ft, fitted2[3]+K1uvt)
print fitted2[4],fitted2[5]
plt.plot (time, ssh, time, double(time,fitted[0], fitted[1], fitted[2], fitted[3]),
time, doublenl(time,fitted2[0],fitted2[1],fitted2[2],fitted2[3],fitted2[4],fitted2[5]))
plt.title(location[stn])
Victoria
M2 corrected 0.289 m and 315 degrees.
M2 better? 0.289 m and 315 degrees.
K1 corrected 0.685 m and 276 degrees.
K1 better? 0.685 m and 276 degrees.
0.0160016867368 -44.4440945659
Out[31]:
<matplotlib.text.Text at 0x4ca4450>
Including an M2 frequency makes no change at the order of mm or degrees.
In [ ]: