Check on a 40 day run with no stratification. This means no OBC temp/salinity, no weather and no rivers. Temp/salinity set to ocean values through whole domain.
In [1]:
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import pandas as pd
import csv
%matplotlib inline
First check that it really is no stratty
In [2]:
path = '/data/nsoontie/MEOPAR/SalishSea/results/stratification/no_strat/2d_grid_T.nc'
gridT=nc.Dataset(path)
In [3]:
temp=gridT.variables['votemper']
sal=gridT.variables['vosaline']
sal.shape
Out[3]:
(20, 40, 898, 398)
In [4]:
fig,axs=plt.subplots(1,2,figsize=(10,10))
ax=axs[0]
mesh=ax.pcolormesh(sal[4,0,:,:])
fig.colorbar(mesh,ax=ax)
ax.set_title('Salinity')
ax=axs[1]
mesh=ax.pcolormesh(temp[4,0,:,:])
fig.colorbar(mesh,ax=ax)
ax.set_title('Temperature')
Out[4]:
<matplotlib.text.Text at 0x7fd82416d350>
Thalweg
In [5]:
dep_d = gridT.variables['deptht']
# Call thalweg
lines = np.loadtxt('/data/nsoontie/MEOPAR/tools/bathymetry/thalweg_working.txt',
delimiter=" ", unpack=False)
lines = lines.astype(int)
ds=np.arange(0,lines.shape[0],1);
XX,ZZ = np.meshgrid(ds,-dep_d[:])
# Salinity along thalweg
salP=sal[-1,:,lines[:,0],lines[:,1]]
salP= np.ma.masked_values(salP,0)
tempP=temp[-1,:,lines[:,0],lines[:,1]]
tempP= np.ma.masked_values(tempP,0)
fig,axs=plt.subplots(2,1,figsize=(20,5))
mesh=axs[0].pcolormesh(XX,ZZ,salP)
axs[0].set_title('salinity')
fig.colorbar(mesh,ax=axs[0])
mesh=axs[1].pcolormesh(XX,ZZ,tempP)
axs[1].set_title('temperature')
fig.colorbar(mesh,ax=axs[1])
Out[5]:
<matplotlib.colorbar.Colorbar instance at 0x7fd80bcf45a8>
Ok, so we had a thermocline... Double check if there is a better way to remove it. Can really see the tidal mixing.
Tides stuff now
In [6]:
# pathname for data - all of the tide runs are stored in this directory
path = '/data/nsoontie/MEOPAR/SalishSea/results/stratification/no_strat/'
#path = '../../myResults/'
#joining the two string together
name = path
print name
/data/nsoontie/MEOPAR/SalishSea/results/stratification/no_strat/
In [7]:
#observations
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 .. ... ... ... ... ... ... ... 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 60 Big Bay 50.36000 125.1300 75.5 14.9 83.3 283.5 61 Chatham Point 50.33000 125.4400 90.3 305.1 65.4 270.5 62 Yorke Island 50.44000 125.9700 117.1 271.8 55.8 260.0 63 Powell River 49.86000 124.5500 100.7 34.3 90.4 286.6 64 Lund 49.98000 124.7600 102.2 35.4 88.9 287.9 65 Nymphe Cove 50.13000 125.3600 61.5 350.4 77.0 279.9 66 Brown Bay 50.16000 125.3700 93.5 315.9 67.9 270.1 67 Maude Island E 50.13000 125.3300 55.6 7.4 81.1 283.9 68 Welsford Island 50.22000 125.1300 99.4 35.1 91.1 286.9 69 Redonda Bay 50.26000 124.9900 97.5 36.7 87.1 287.4 70 Channel Islands 50.31000 124.7500 102.6 35.9 89.9 288.0 71 Turnback Point 50.42000 125.1200 102.0 37.0 91.7 287.6 72 Orford Bay 50.59000 124.8600 101.5 37.2 90.3 288.1 73 Waddington Harbour 50.87000 124.8700 103.4 38.0 89.2 288.2 74 Shoal Bay 50.46000 125.3600 89.9 307.5 66.6 269.6 75 Kelsey Bay 50.39000 125.9600 117.0 276.3 57.7 261.4 76 Tacoma 47.26670 122.4133 113.9 11.8 83.8 277.9 [77 rows x 7 columns]
In [10]:
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']
In [27]:
stn_lon=np.zeros(numsta); stn_lat=np.zeros(numsta)
#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]
stn_lon[stn]=lon
stn_lat[stn]=lat
In [12]:
#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 April 21 2003
# data for phase output from bdytides.F90; found in ocean.output
K1ft = 1.065505
K1uvt = 111.481741
M2ft = 0.982328
M2uvt = 250.506179
O1ft = 1.105495
O1uvt = 142.040782
S2ft = 1.0
S2uvt = 0.0
P1ft = 1.0
P1uvt =241.335269
N2ft = 0.982328
N2uvt =205.684028
Q1ft = 1.105495
Q1uvt =97.218631
K2ft = 1.159036
K2uvt = 42.361669
In [13]:
# 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.))
In [14]:
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)
PortRenfrew Victoria PointAtkinson CampbellRiver
In [21]:
#allocate space for our arrays
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
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
if fitted[2] < 0:
fitted[2] = -fitted[2]
fitted[3] = fitted[3]+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)
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]
350.472573195
In [22]:
def mean(diff):
return np.mean(abs(diff))
In [23]:
def rms(diff):
return np.sqrt(np.mean(diff**2))
In [24]:
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
In [25]:
#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))
In [26]:
outfile = 'tides_nostrat.csv'
split1=8; split2=20
with open(outfile, 'wb') as csvfile:
writer = csv.writer(csvfile, delimiter=',')
writer.writerow([
'Station Number', 'Station Name', 'Longitude', 'Latitude',
'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( [str(stn +1),stations_obs[stn], str(stn_lon[stn]), str(stn_lat[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])])
In [ ]: