#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)


#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')

# 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)

#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)