Compare harmonics from currents between measured and modelled. Measured harmonics are taken from Table 3 of Foreman et al (1995)
#Preamble
%matplotlib inline
from salishsea_tools import tidetools
#Load in text file containing measured current constituents taken from Table 3 of Foreman et al (1995).
#latitude and longitude taken from http://www.charts.gc.ca/copyright-droitdauteur/documents/currents-courants-eng.asp
import pandas
import numpy as np
from math import radians, cos, sin
measured = pandas.read_csv('/ocean/klesouef/meopar/tools/compare_tides/obs_tidal_current_const_Foreman95.csv'\
,delimiter=';')#,index_col='Site')
#convert M2 velocity to U and V component, and from cm/s to m/s
M2_amp_U = measured.M2amp*np.cos(np.radians(measured.M2dir))/100
M2_amp_V = measured.M2amp*np.sin(np.radians(measured.M2dir))/100
#find the model co-ordinates of each of the measured locations
model_x1 = []
model_y1 = []
for i in np.arange(0,len(measured.Lon)):
x1, y1 = tidetools.find_closest_model_point(-measured.Lon[i],measured.Lat[i])
model_x1.append(x1)
model_y1.append(y1)
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-3-aeb5aab95bb9> in <module>() 4 5 for i in np.arange(0,len(measured.Lon)): ----> 6 x1, y1 = tidetools.find_closest_model_point(-measured.Lon[i],measured.Lat[i]) 7 model_x1.append(x1) 8 model_y1.append(y1) TypeError: find_closest_model_point() takes exactly 5 arguments (2 given)
#Get modelled current harmonics
import netCDF4 as NC
harmU = NC.Dataset('/ocean/klesouef/meopar/tools/NetCDF_Plot/WC3_Harmonics_gridU_TIDE2D.nc','r')
harmV = NC.Dataset('/ocean/klesouef/meopar/tools/NetCDF_Plot/WC3_Harmonics_gridV_TIDE2D.nc','r')
X_U = harmU.variables['nav_lon']
Y_U = harmU.variables['nav_lat']
#Get amplitude and phase
mod_U_M2_amp = harmU.variables['M2_amp'][0,:,:]
mod_U_M2_pha = harmU.variables['M2_pha'][0,:,:]
mod_V_M2_amp = harmV.variables['M2_amp'][0,:,:]
mod_V_M2_pha = harmV.variables['M2_pha'][0,:,:]
#print a text file showing the modelled and measured components (similar to Table 3 of Foreman et al 1995)
import csv
loc = "/ocean/klesouef/meopar/tools/compare_tides"
with open('current_harm_comps.csv', 'wb') as csvfile:
writer = csv.writer(csvfile, delimiter=',')
writer.writerow(['Station Number','Station Name','Longitude','Latitude','Model U M2 Amp','Model U M2 Phase',\
'Model V M2 Amp','Model V M2 Phase','Measured U M2 Amp','Measured U M2 Phase',\
'Measured V M2 Amp','Measured V M2 Phase'])
for t in np.arange(0,len(measured.Lat)):
a = model_x1[t]
b = model_y1[t]
print(a,b)
writer.writerow([str(t+1),measured.Site[t],-measured.Lon[t],measured.Lat[t],mod_U_M2_amp[a,b],mod_U_M2_pha[a,b],\
mod_V_M2_amp[a,b],mod_V_M2_pha[a,b],M2_amp_U[t],'??',M2_amp_V[t],'??'])
print('Results saved here: '+loc+'/current_harm_diffs.csv')
(76, 150) (89, 159) (181, 262) (230, 240) (261, 235) (266, 222) Results saved here: /ocean/klesouef/meopar/tools/compare_tides/current_harm_diffs.csv
# get the bathy to check if this grid point is on land ('masked')
grid = NC.Dataset('/ocean/klesouef/meopar/nemo-forcing/grid/SubDom_bathy_meter_NOBCchancomp.nc','r')
bathy = grid.variables['Bathymetry'][:,:]
# define a function to plot the measured and model locations
def plot_mod_meas_on_bathy(t):
modlon = X_U[model_x1[t],model_y1[t]]
modlat = Y_U[model_x1[t],model_y1[t]]
measlon = -measured.Lon[t]
measlat = measured.Lat[t]
plt.contourf(X_U,Y_U,bathy)
plt.colorbar()
plt.title('Domain of model (depths in m)')
hold = True
plt.plot(modlon,modlat,'g.',markersize=10,label='model')
plt.plot(measlon,measlat,'m.',markersize=10,label='measured')
plt.xlim([modlon-0.2,modlon+0.2])
plt.ylim([modlat-0.2,modlat+0.2])
plt.legend(numpoints=1)
Some of the measured locations don't have any model points near them
#Active Pass is not well resolved in the bathymetry
plot_mod_meas_on_bathy(2)
#Gabriola Passage is not well resolved in the bathymetry
plot_mod_meas_on_bathy(4)
#Dodd Narrows is not well resolved in the bathymetry
plot_mod_meas_on_bathy(5)