# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 2
#     language: python
#     name: python2
# ---

# This notebook will be used to look at and understand the model output currents through Haro Strait with the use of ellipses

# +
import os
import datetime
import matplotlib.pylab as plt
from matplotlib.patches import Ellipse
import numpy as np
from IPython.display import display, Math, Latex
import csv
import pandas as pd
import subprocess as sp
from subprocess import Popen, PIPE
import seaborn as sns

import netCDF4 as nc
from salishsea_tools import viz_tools

# %matplotlib inline
# -

font = {'family' : 'Arial',
    'weight' : 'normal',
    'size'   : 20}
axisfont = {'family' : 'Arial',
    'weight' : 'light',
    'size'   : 16}
sns.set_style('darkgrid')

# #Which grid points do we want to focus on?

thalweg = np.loadtxt('/data/dlatorne/MEOPAR/tools/bathymetry/thalweg_working.txt', dtype=int, unpack=True)
grid = '/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc'
fB = nc.Dataset(grid)
bathy = fB.variables['Bathymetry'][:]
bathy = np.ma.masked_values(bathy,0)
lons = fB.variables['nav_lon'][:]
lats = fB.variables['nav_lat'][:]

#Haro Strait
iH = 233
jH = 312

cmap = plt.get_cmap('winter_r')
cmap.set_bad('burlywood')
fig, ax = plt.subplots(1,1)
plt.pcolormesh(bathy, cmap=cmap)
plt.plot(iH,jH,'ro')
plt.plot(thalweg[1,480:520:4],thalweg[0,480:520:4], 'ko')
plt.axis([210,250,280,360])
ax.set_title('Haro Strait', **font)
ax.set_ylabel('y index', **axisfont)
ax.set_xlabel('x index', **axisfont)

lats[jH,iH], lons[jH, iH]

# +
fig, ax = plt.subplots(1,1)
viz_tools.set_aspect(ax) 
ax.contourf(lons,lats,bathy.data,cmap=plt.get_cmap('winter_r'))
ax.plot(lons[jH, iH], lats[jH, iH],'ro')
ax.plot(lons[thalweg[0,480:530:4],thalweg[1,480:530:4]], lats[thalweg[0,480:530:4],thalweg[1,480:530:4]], 'mo')
viz_tools.plot_land_mask(ax, fB, coords='map')
viz_tools.plot_coastline(ax, fB, coords='map', isobath=5)

ax.set_xlim([lons[jH, iH]-0.12, lons[jH, iH]+0.12])
ax.set_ylim([lats[jH, iH]-0.1, lats[jH, iH]+0.18])
ax.set_ylabel('Latitude (degrees N)', **axisfont)
ax.set_xlabel('Longitude (degrees W)', **axisfont)
ax.set_title('Haro Strait', **font)
ax.tick_params(axis='both', which='major', labelsize=14)


# +
to = datetime.datetime(2014,11,26)
tf = datetime.datetime(2015, 7, 13)

t_o = to.strftime('%d%b%y').lower()
t_f = tf.strftime('%d%b%y').lower()

loc = np.arange(480,520,4)
count = np.arange(0,len(loc))
parh = np.zeros((len(loc),39,8))
parh_av = np.zeros((len(loc),8))
cols=np.arange(0,8)

i = np.zeros((16))
j = np.zeros((16))


for k,l in zip(loc,count):
    i = thalweg[1,loc]
    j = thalweg[0,loc]
    runname1 = '{}_{}_HaroS_{}'.format(t_o, t_f, k)
    runname2 = '{}_{}_HaroS_{}_depav(20-400)'.format(t_o, t_f, k)
    df1 = pd.read_csv('/ocean/mdunn/MEOPAR/analysis/Muriel/TidalEllipseData/HaroStrait/'+runname1+'.csv', usecols=cols) 
    parh[l,:,:] = df1.as_matrix()
    df2 = pd.read_csv('/ocean/mdunn/MEOPAR/analysis/Muriel/TidalEllipseData/HaroStrait/'+runname2+'.csv', index_col=False)
    parh_av[l,:] = df2.as_matrix()
# -

new = np.zeros((parh.ndim-1,len(parh)))
new.shape

# ##Model Ellipses with Depth

# +
# Set up the figure and axes
count = np.arange(0,len(loc))
loc = np.arange(480,520,4)
fig, (axl, axcb) = plt.subplots(1, 2, figsize=(20, 10))
land_colour = 'burlywood'
axl.set_axis_bgcolor(land_colour)
axl.set_position((0.125, 0.125, 0.6, 0.775))
axcb.set_position((0.73, 0.125, 0.02, 0.775))

smin, smax, dels = 29, 35, 0.5
cmap = plt.get_cmap('Blues')
cmap.set_bad(land_colour)
tracers = nc.Dataset('/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/18sep27sep/SalishSea_1d_20030918_20030927_grid_T.nc')
sal = tracers.variables['vosaline']
npsal = sal[:]
zlevels = tracers.variables['deptht']
sal_0 = npsal[-1, :, thalweg[0], thalweg[1]]
sal_tzyx = np.ma.masked_values(sal_0, 0)
x, z = np.meshgrid(np.arange(thalweg.shape[1]), zlevels)
mesh = axl.pcolormesh(x, z, sal_tzyx.T, cmap=cmap, vmin=smin, vmax=smax)
cbar = plt.colorbar(mesh, cax=axcb)
cbar.set_label('Practical Salinity')

scale = 10
for c, ind in zip(count, loc):
    for dep, value in zip(np.arange(0,39), parh[0,:,0]):
        if parh[c,dep,2] > 0:
            thec = 'b'
        else:
            thec = 'r'
        ellsc = Ellipse(xy=(ind, value), width=scale*parh[c,dep,1], 
                            height=scale*parh[c,dep,2], angle=parh[c, dep,3],
                            color=thec)
        axl.add_artist(ellsc)
        ellsc.set_facecolor(thec)

axl.set_ylim(-2,400)
axl.invert_yaxis()
axl.set_xlim(475,525)
axl.set_title('Tidal ellipses in Haro Strait')
axl.set_ylabel('Depth (m)')
axl.set_xlabel('Distance along the thalweg')
# -

# They are mostly going northwards... larger values at the surface and at depth and smaller in the middle.

# ##Load Foreman values

# +
i = thalweg[1,loc]
j = thalweg[0,loc]

lon = lons[j,i]
lat = lats[j,i]
print lat, lon.shape
# -

parhf = np.zeros((10,8))
os.chdir(r'/ocean/mdunn/Tides/')
for k in count:
    p = Popen(["matlab","-nodesktop","-nodisplay","-r", "findparam(%s,%s)" % (lon[k], lat[k])], stdout=PIPE, stderr=PIPE)
    output, err = p.communicate()
    parhf[k] = np.loadtxt('output.txt')

# ##Load current meter locations parameters

csv1 = pd.read_csv('/ocean/mdunn/MEOPAR/analysis/Muriel/param_mat.csv', index_col=False)
param_mat = csv1.as_matrix()
csv2 = pd.read_csv('/ocean/mdunn/MEOPAR/analysis/Muriel/params.csv', index_col=False)
params = csv2.as_matrix()

# ##Map of ellipses

# +
phi=29
fig = plt.figure(figsize=(20,15))
k = np.zeros((898,398)); m = np.zeros((898,398))
ax = fig.add_subplot(111)
viz_tools.set_aspect(ax) 
ex = 80
scale = 5
imin = min(i)
imax = max(i)

jmin = min(j)
jmax = max(j)

for q in np.arange(jmin-ex,jmax+ex):
    for l in np.arange(imin-ex,imax+ex):
        k[q,l] = q*np.cos(phi*np.pi/180.)+l*np.sin(phi*np.pi/180.)
        m[q,l] = -q*np.sin(phi*np.pi/180.)+l*np.cos(phi*np.pi/180.) 

for index, x, y in zip(count, i ,j):
    if parh_av[index,1] > 0:
        thec = 'b'
    else:
        thec = 'r'
    ellsc = Ellipse(xy=(m[y,x],k[y,x]), width=scale*parh_av[index,0], 
                        height=scale*parh_av[index,1], angle=parh_av[index,2],
                        color=thec)
    ax.add_artist(ellsc)
    ellsc.set_facecolor(thec)
    
    if parhf[index,1] > 0:
        thec = [0,0,0.6, 0.3]
    else:
        thec = [0.6,0,0, 0.3]
    ellsc = Ellipse(xy=(m[y,x],k[y,x]), width=scale*parhf[index,0],
                        height=scale*parhf[index,1], angle=parhf[index,2],
                        color=thec)
    ax.add_artist(ellsc)
    ellsc.set_facecolor(thec)
    
#Plot the Haro Strait current meter location    
ax.plot(m[jH,iH], k[jH,iH],'g^', markersize=12)   

#Values from EllipsesComparisonswithForemanModel notebook

if params[2,1] > 0:
    thec = 'b'
else:
    thec = 'r'
ellsc = Ellipse(xy=(m[jH,iH],k[jH,iH]), width=scale*params[2,0], 
                    height=scale*params[2,1], angle=params[2,2],
                    color=thec)
ax.add_artist(ellsc)
ellsc.set_facecolor(thec)

if param_mat[2,1] > 0:
    thec = [0,0,0.6, 0.3]
else:
    thec = [0.6,0,0, 0.3]
ellsc = Ellipse(xy=(m[jH,iH],k[jH,iH]), width=scale*param_mat[2,0],
                    height=scale*param_mat[2,1], angle=param_mat[2,2],
                    color=thec)
ax.add_artist(ellsc)
ellsc.set_facecolor(thec)

# y-axis in k, but labelled in latitude
ax.set_ylim(375,425)
slope = (lats[jmax,imax]-lats[jmin,imin])/(k[jmax,imax]-k[jmin,imin])
mylist = (k[jmax,imax]+(np.arange(48.335, 48.58, 0.05)-
                           lats[jmin,imin])/slope).tolist()
labels = ['48.5','48.55', '48.6', '48.65', '48.7', '48.75']
ax.set_yticks(mylist)
ax.set_yticklabels(labels)
ax.set_ylabel('Latitude (degrees N)')
# x-axis in m, but labelled in longitude
ax.set_xlim(30,80)
slope = (lons[jmax,imax]-lons[jmin,imin])/(m[jmax,imax]-m[jmin,imin])
mylist = (m[jmax,imax]+(np.arange(-123.4,-123,0.1)-
                        lons[jmin,imin])/slope).tolist()
labels = ['123.4','123.3', '123.2', '123.1', '123']
ax.set_xticks(mylist)
ax.set_xticklabels(labels)
ax.set_xlabel('Longitude (degrees W)')

# land, and 5 m contour
contour_interval = [-0.01, 0.01]
ax.contourf(m[jmin-ex:jmax+ex,imin-ex:imax+ex],k[jmin-ex:jmax+ex,imin-ex:imax+ex],
        bathy.data[jmin-ex:jmax+ex,imin-ex:imax+ex],contour_interval,colors='black')
ax.contour(m[jmin-ex:jmax+ex,imin-ex:imax+ex],k[jmin-ex:jmax+ex,imin-ex:imax+ex],
            bathy.data[jmin-ex:jmax+ex,imin-ex:imax+ex],[5],colors='black')
ax.set_title('Depth averaged M2 tidal ellipses through Haro Strait')

print "red is clockwise"
# -

# ##Line graph with observations

parameter = ['major-axis', 'minor-axis', 'inclination']
unit = ['m/s' ,'m/s','deg']
obs = np.ones((10))*0.513
obspha = [np.ones((10))*321,np.ones((10))*221]
Hlat = np.ones((2))*48.537579
y = [0.4,1.2]

# +
fig, (ax0, ax1, ax2, ax3, ax4)= plt.subplots(5,1,figsize=(20,20), sharex=True)

for ax , k, par, u in zip([ax0, ax1, ax2], [0,1,2], parameter, unit):
    ax.plot(lat[:],parh_av[:,k], '-', label='Salish Sea Model')
    ax.plot(lat[:], parhf[:,k], '-', label='Foreman Model')
    ax.set_title('Comparing the {}.'.format(par), **font)
    ax.set_ylabel('{} ({})'.format(par,u),**axisfont)
    ax.tick_params(axis='both', which='major', labelsize=14)
ax0.plot(Hlat, y, 'r--')
ax0.set_ylim([0.47,1.1])
ax0.plot(lat[:], obs, label='Observations by CM')
ax0.plot(48.537579, params[2,0], 'bo', label = 'Salish Sea Model at CM')
ax0.plot(48.537579, param_mat[2,0], 'go', label = 'Foreman Model at CM')
ax0.legend(loc=0, fontsize=16)

for ax, ind, k in zip([ax3,ax4], [0,1], [3,7]):
    ax.plot(lat[:],parh_av[:,k], '-', label='Salish Sea Model')
    ax.plot(lat[:], parhf[:,k], '-', label='Foreman Model')
    ax.plot(lat[:], obspha[ind], label='Observations by CM')
    ax.plot(Hlat[0], params[2,k], 'bo', label = 'Salish Sea Model at CM')
    ax.plot(Hlat[0], param_mat[2,k], 'go', label = 'Foreman Model at CM')
    ax.plot(Hlat, [0,400], 'r--')
    ax.set_ylabel('Phase (deg)', **axisfont)
    ax.set_ylim([0,370])
    ax.tick_params(axis='both', which='major', labelsize=14)
ax4.set_xlabel('Latitude', **axisfont)
ax3.set_title('Comparing the M2 Phase', **font)
ax4.set_title('Comparing the K1 Phase', **font)
# -





StarcoderdataJupyterScriptsDedupFiltered