This notebook generates tidal comparisons for CMOS.
# 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
import pandas as pd
import csv
# grid
bathy, X, Y = tidetools.get_SS_bathy_data()
# pathname for data
#This run has bf=1e-3, K1 phase and amplitude corrected and M2 phase corrected
name = '/data/nsoontie/MEOPAR/SalishSea/results/tides/tide_M2phase/'
#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.
# 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
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).
filename = '/data/nsoontie/MEOPAR/analysis/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]
Here is a list of the stations for which we have model results.
stations = ['BrownBay', 'CampbellRiver', 'ChathamPoint', 'CloverPoint',
'Esquimalt', 'FinnertyCove', 'FulfordHarbour', 'GibsonsLanding',
'HalfmoonBay', 'IrvinesLanding', 'KelseyBay', 'Lund',
'MaudeIslandE', 'NympheCove', 'PatosIsland','PedderBay',
'PointAtkinson','PointGrey','PortRenfrew',
'PowellRiver', 'Sandheads','SeymourNarrows','SheringhamPoint',
'Tsawwassen','TumboChannel','TwinIslets','Victoria',
'WhalerBay','YorkeIsland']
numsta=len(stations)
#again with spaces because the text file likes that
stations_obs = ['Brown Bay', 'Campbell River', 'Chatham Point', 'Clover Point',
'Esquimalt', 'Finnerty Cove', 'Fulford Harbour', 'Gibsons Landing',
'Halfmoon Bay','Irvines Landing', 'Kelsey Bay', 'Lund',
'Maude Island E', 'Nymphe Cove', 'Patos Island','Pedder Bay',
'Point Atkinson','Point Grey','Port Renfrew',
'Powell River', 'Sandheads','Seymour Narrows','Sheringham Point',
'Tsawwassen','Tumbo Channel','Twin Islets','Victoria',
'Whaler Bay','Yorke Island']
#Removed Port Townsend because it was recorded twice in the file. Which to use?
M2_amp=[]
M2_pha=[]
K1_amp=[]
K1_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)
ts = 240
for stn in range(numsta):
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]
fitted, cov = curve_fit(double,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)
#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]
BrownBay CampbellRiver ChathamPoint CloverPoint Esquimalt FinnertyCove FulfordHarbour GibsonsLanding HalfmoonBay IrvinesLanding KelseyBay Lund MaudeIslandE NympheCove PatosIsland PedderBay PointAtkinson PointGrey PortRenfrew PowellRiver Sandheads SeymourNarrows SheringhamPoint Tsawwassen TumboChannel TwinIslets Victoria WhalerBay YorkeIsland
print M2_amp_obs
print K1_amp_obs
[ 0.935 0.825 0.903 0.403 0.367 0.447 0.582 0.947 0.964 0.988 1.17 1.022 0.556 0.615 0.68 0.342 0.918 0.945 0.708 1.007 0.869 0.946 0.487 0.811 0.726 1.013 0.373 0.834 1.171] [ 0.679 0.846 0.654 0.642 0.643 0.708 0.753 0.872 0.88 0.88 0.577 0.889 0.811 0.77 0.79 0.627 0.862 0.906 0.453 0.904 0.837 0.692 0.548 0.834 0.811 0.904 0.627 0.847 0.558]
#Plotting - M2
font = {
'size' : 18}
plt.rc('font', **font)
fig=tidetools.plot_scatter_pha_amp(M2_amp,M2_amp_obs,M2_pha,M2_pha_obs,'M2',figsize=(14,6))
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)
[<matplotlib.lines.Line2D at 0x4f262d0>]
#Plotting - K1
font = {
'size' : 16}
plt.rc('font', **font)
fig=tidetools.plot_scatter_pha_amp(K1_amp,K1_amp_obs,K1_pha,K1_pha_obs,'K1',figsize=(14,6))
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)
[<matplotlib.lines.Line2D at 0x38c4c50>]
# Port Townsend
stn='PortTownsend'
fT1 = NC.Dataset(name+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
M2_ampPT=fitted[0]*M2ft
pha = fitted[1]+M2uvt
if pha > 360:
pha=pha-360
M2_phaPT=pha
K1_ampPT=fitted[2]*K1ft
pha= fitted[3]+K1uvt
if pha > 360:
pha=pha-360
K1_phaPT=pha
print M2_ampPT,M2_phaPT,K1_ampPT, K1_phaPT
0.537417179642 343.912392102 0.750546673222 270.711264264
At Port Townsend, M2 is pretty off from both Foreman et al 1995 and 2004. (Closer to Foreman et al 1995 )
def rms(diff):
return np.sqrt((diff**2).mean())
diff =M2_amp-M2_amp_obs
print 'M2 amp mean error', np.mean(abs(diff))
print 'M2 amp rms error', rms(diff)
diff =M2_pha-M2_pha_obs
diff = diff +360*(diff<-180) - 360*(diff>180)
print 'M2 pha mean error', np.mean(abs(diff))
print 'M2 pha rms error', rms(diff)
M2 amp mean error 0.0712630845562 M2 amp rms error 0.125563550289 M2 pha mean error 6.8524786232 M2 pha rms error 13.1054555665
diff =K1_amp-K1_amp_obs
print 'K1 amp mean error', np.mean(abs(diff))
print 'K1 amp rms error', rms(diff)
diff =K1_pha-K1_pha_obs
diff = diff +360*(diff<-180) - 360*(diff>180)
print 'K1 pha mean error', np.mean(abs(diff))
print 'K1 pha rms error', rms(diff)
K1 amp mean error 0.02552682652 K1 amp rms error 0.0447435264431 K1 pha mean error 2.70244129003 K1 pha rms error 4.52104797531
diff =M2_amp-M2_amp_obs
print diff
[-0.12226278 -0.04050405 -0.01621579 -0.17383465 -0.15478292 -0.08726851 -0.09313465 -0.02236256 0.00449838 -0.00519835 0.04330538 0.00559139 -0.06427212 -0.13347857 0.0309048 -0.14615282 0.00803197 -0.01661162 -0.02174928 -0.00504801 0.01896887 -0.54028503 -0.05359559 0.00642456 0.0449969 0.00783555 -0.15361937 0.00809151 0.03760346]
bad station for M2 Amp living in index 21
print stations[21]
SeymourNarrows
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(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
D_M2= complex_diff(np.array(M2_amp_obs),np.array(M2_pha_obs), np.array(M2_amp),np.array(M2_pha))
print 'M2 Complex diff mean:', np.mean(D_M2)
D_K1= complex_diff(np.array(K1_amp_obs),np.array(K1_pha_obs), np.array(K1_amp),np.array(K1_pha))
print 'K1 Complex diff mean:', np.mean(D_K1)
M2 Complex diff mean: 0.115614209695 K1 Complex diff mean: 0.043362009465