#!/usr/bin/env python
# coding: utf-8

# # Check River Climatology

# In[64]:


from __future__ import division

import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np

from salishsea_tools import nc_tools

get_ipython().run_line_magic('matplotlib', 'inline')


# In[65]:


def find_points(flow):
    for i in range(390,435):
        for j in range(280,398):
            if flow[i,j] > 0:
                print i,j, lat[i,j], lon[i,j], flow[i,j]


# In[66]:


grid = nc.Dataset('/ocean/sallen/allen/research/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')
lat = grid.variables['nav_lat'][:,:]
lon = grid.variables['nav_lon'][:,:]
depth = grid.variables['Bathymetry'][:]


# In[67]:


river = nc.Dataset('/ocean/sallen/allen/research/MEOPAR/NEMO-forcing/rivers/rivers_month.nc')
print river.variables['rorunoff'][:].shape


# In[68]:


print 'Apr'
find_points(river.variables['rorunoff'][4,:,:])
print 'May'
find_points(river.variables['rorunoff'][5,:,:])


# In[69]:


ik = 425; jk = 302; d = 6
fig, ax = plt.subplots(1,1,figsize=(15,7.5))
imin = 390; imax = 435; jmin = 280; jmax = 398
cmap = plt.get_cmap('winter_r')
cmap.set_bad('burlywood')
mesh = ax.pcolormesh(depth[imin:imax,jmin:jmax], vmax = 10., cmap=cmap)
ax.set_xlim((0,110))
ax.set_xlabel('Grid Points')
ax.set_ylabel('Grid Points')
ax.text(40, 28, "Short Fraser River", fontsize=14)
cbar=fig.colorbar(mesh)
cbar.set_label('Depth (m)')
ax.plot(np.array((324,324,334,334,334,318))-jmin+0.5,np.array((411,412,414,415,416,434))-imin+0.5,'ko');
ax.plot(np.array((322-1,318-2))-jmin+0.5,np.array((412-1,414-1))-imin+0.5,'rs');


# In[71]:


data = nc.Dataset('/ocean/sallen/allen/research/MEOPAR/myResults/NEMO36_Tides/base_run/E1/SalishSea_1d_20030421_20030424_grid_T.nc')
nc_tools.show_variables(data)


# In[72]:


S = data.variables['vosaline'][:]
print S.shape


# In[73]:


t=3; d = 0
imin = 390; imax = 435; jmin = 280; jmax = 398
lS = np.log10(S[t,d,imin:imax,jmin:jmax])
fig, ax = plt.subplots(1,1,figsize=(15,7.5))
cmap = plt.get_cmap('spectral')
cmap.set_bad('burlywood')
mesh = ax.pcolormesh(lS, vmax=1.2, vmin=-40, cmap=cmap)
ax.set_xlim((0,110))
ax.set_xlabel('Grid Points')
ax.set_ylabel('Grid Points')
ax.text(40, 28, "Short Fraser River", fontsize=14)
cbar=fig.colorbar(mesh)
cbar.set_label('Salinity (practical)')
ax.plot(np.array((324,324,334,334,334,318))-jmin+0.5,np.array((411,412,414,415,416,434))-imin+0.5,'wo');
ax.plot(np.array((322-1,318-2))-jmin+0.5,np.array((412-1,414-1))-imin+0.5,'rs');


# In[74]:


datau = nc.Dataset('/ocean/sallen/allen/research/MEOPAR/myResults/NEMO36_Tides/base_run/E1/SalishSea_1d_20030421_20030424_grid_U.nc')
u = datau.variables['vozocrtx'][:]


# In[75]:


t = 3; d = 0
imin = 390; imax = 435; jmin = 280; jmax = 398
fig, ax = plt.subplots(1,1,figsize=(15,7.5))
cmap = plt.get_cmap('spectral')
cmap.set_bad('burlywood')
mesh = ax.pcolormesh(u[t,d,imin:imax,jmin:jmax], cmap=cmap)
ax.set_xlim((0,110))
ax.set_xlabel('Grid Points')
ax.set_ylabel('Grid Points')
ax.text(40, 28, "Short Fraser River", fontsize=14)
cbar=fig.colorbar(mesh)
cbar.set_label('Cross-page velocity')
ax.plot(np.array((324,324,334,334,334,318))-jmin+0.5,np.array((411,412,414,415,416,434))-imin+0.5,'wo');
ax.plot(np.array((322-1,318-2))-jmin+0.5,np.array((412-1,414-1))-imin+0.5,'rs');


# In[76]:


datav = nc.Dataset('/ocean/sallen/allen/research/MEOPAR/myResults/NEMO36_Tides/base_run/E1/SalishSea_1d_20030421_20030424_grid_V.nc')
v = datav.variables['vomecrty'][:]


# In[77]:


t = 3; d = 0
imin = 390; imax = 435; jmin = 280; jmax = 398
fig, ax = plt.subplots(1,1,figsize=(15,7.5))
cmap = plt.get_cmap('spectral')
cmap.set_bad('burlywood')
mesh = ax.pcolormesh(v[t,d,imin:imax,jmin:jmax], cmap=cmap)
ax.set_xlim((0,110))
ax.set_xlabel('Grid Points')
ax.set_ylabel('Grid Points')
ax.text(40, 28, "Short Fraser River", fontsize=14)
cbar=fig.colorbar(mesh)
cbar.set_label('Up-page velocity')
ax.plot(np.array((324,324,334,334,334,318))-jmin+0.5,np.array((411,412,414,415,416,434))-imin+0.5,'wo');
ax.plot(np.array((322-1,318-2))-jmin+0.5,np.array((412-1,414-1))-imin+0.5,'rs');


# In[78]:


ssh =  data.variables['sossheig'][:]


# In[79]:


imin = 390; imax = 435; jmin = 280; jmax = 398
fig, ax = plt.subplots(1,1,figsize=(15,7.5))
cmap = plt.get_cmap('spectral')
cmap.set_bad('burlywood')
mesh = ax.pcolormesh(ssh[t,imin:imax,jmin:jmax], cmap=cmap)
ax.set_xlim((0,110))
ax.set_xlabel('Grid Points')
ax.set_ylabel('Grid Points')
ax.text(40, 28, "Short Fraser River", fontsize=14)
cbar=fig.colorbar(mesh)
cbar.set_label('Sea Surface Height')
ax.plot(np.array((324,324,334,334,334,318))-jmin+0.5,np.array((411,412,414,415,416,434))-imin+0.5,'wo');
ax.plot(np.array((322-1,318-2))-jmin+0.5,np.array((412-1,414-1))-imin+0.5,'rs');


# In[80]:


# point of interest is 413, 317
print S[0,0,413,317], S[0,1,413,317]


# In[81]:


print u[0,0,413,317], S[0,0,413,318]


# In[82]:


print u[0,0,413,316], S[0,0,413,316]


# In[83]:


print v[0,0,413,317], S[0,0,414,317]


# In[84]:


print v[0,0,412,317], S[0,0,412,317]


# In[85]:


print -u[0,0,413,317]*S[0,0,413,318]+u[0,0,413,316]*S[0,0,413,316]
print -v[0,0,413,317]*S[0,0,414,317]+v[0,0,412,317]*S[0,0,412,317]


# In[ ]: