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

# Notebook to explore Live Ocean results from Parker MacCready and group (http://faculty.washington.edu/pmacc/LO/LiveOcean.html). Goal is to determine a region for extracting boundary conditions.

# In[1]:


import netCDF4 as nc
import matplotlib.pyplot as plt
import UBC_subdomain

from salishsea_tools import viz_tools

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


# In[2]:


f = nc.Dataset('/ocean/nsoontie/MEOPAR/LiveOcean/20160608/ocean_his_0025.nc')


# In[3]:


f.variables


# Lots of output! I am interested in:
# * salt
# * temp (which is described as potential temperature)
# * lat_rho, lon_rho, mask_rho
# * pn, pm (1/DY, 1/DX)
# * ocean_time
# * s_rho, hc, Cs_r, Vtransform (vertical grid stuff: Parker sent some code for translating into a depth)
# * zeta (ssh)
# * h (bathymetry at rho points)
# 
# Anything else?

# ### Timing variables

# In[4]:


dt = f.variables['dt']
ntimes = f.variables['ntimes']
ntimes[:]*dt[:]/86400


# It appears as if this file was part of a three day simulation.

# In[5]:


ocean_time = f.variables['ocean_time']
ocean_time[:]


# In[6]:


ocean_time


# In[7]:


time = nc.num2date(ocean_time[:], ocean_time.units)
time


# It is fairly easy to grab the snap shot time using ocean_time and the nc.num2date method.
# 
# ### Surface plots of salinity

# In[8]:


sal = f.variables['salt']
lon_rho = f.variables['lon_rho']
lat_rho = f.variables['lat_rho']
s_rho = f.variables['s_rho']


# In[9]:


sal


# In[10]:


b = nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')
X = b.variables['nav_lon']
Y = b.variables['nav_lat']


# In[11]:


fig,ax=plt.subplots(1,1)
mesh=ax.pcolormesh(lon_rho, lat_rho, sal[0,0,:,:])
cbar=plt.colorbar(mesh,ax=ax)
cbar.set_label('salinity')
viz_tools.plot_coastline(ax,b,coords='map')
iss = 55
iee = 80
jss = 295
jee = 325
ax.plot([lon_rho[jss,iss], lon_rho[jss,iee]],
         [lat_rho[jss,iss], lat_rho[jss,iee]], '-k')
ax.plot([lon_rho[jss,iss], lon_rho[jee,iss]],
         [lat_rho[jss,iss], lat_rho[jee,iss]], '-k')
ax.plot([lon_rho[jee,iss], lon_rho[jee,iee]],
         [lat_rho[jee,iss], lat_rho[jee,iee]], '-k')
ax.plot([lon_rho[jee,iee], lon_rho[jss,iee]],
         [lat_rho[jee,iee], lat_rho[jss,iee]], '-k')


# In[12]:


fig,ax=plt.subplots(1,1)
mesh=ax.pcolormesh(lon_rho[jss:jee,iss:iee],
               lat_rho[jss:jee,iss:iee],
               sal[0,0,jss:jee,iss:iee])
cbar=plt.colorbar(mesh,ax=ax)
cbar.set_label('salinity')
viz_tools.plot_coastline(ax,b,coords='map')
ax.set_xlim([-125.1,-124.4])
ax.set_ylim([48.3,48.9])
ax.plot(X[:,10], Y[:,10],'g-')


# ### Grid and staggering

# In[13]:


lon_u = f.variables['lon_u']
lat_u = f.variables['lat_u']
lon_v = f.variables['lon_v']
lat_v = f.variables['lat_v']


# In[14]:


print('Bottom left')
print('U', lon_u[0,0], lat_u[0,0])
print('V', lon_v[0,0], lat_v[0,0])
print('sal', lon_rho[0,0], lat_rho[0,0])


# In[15]:


print('Top left')
print('U', lon_u[-1,0], lat_u[-1,0])
print('V', lon_v[-1,0], lat_v[-1,0])
print('sal', lon_rho[-1,0], lat_rho[-1,0])


# In[16]:


print('Bottom right')
print('U', lon_u[0,-1], lat_u[0,-1])
print('V', lon_v[0,-1], lat_v[0,-1])
print('sal', lon_rho[0,-1], lat_rho[0,-1])


# In[17]:


print('Top right')
print('U', lon_u[-1,-1], lat_u[-1,-1])
print('V', lon_v[-1,-1], lat_v[-1,-1])
print('sal', lon_rho[-1,-1], lat_rho[-1,-1])


# ### Bathymetry

# In[18]:


h = f.variables['h']
h


# In[19]:


fig,ax=plt.subplots(1,1)
mesh=ax.pcolormesh(lon_rho, lat_rho, h,vmin=0,vmax=300)
cbar=plt.colorbar(mesh,ax=ax)
cbar.set_label('bathy [m]')
viz_tools.plot_coastline(ax,b,coords='map')
ax.plot([lon_rho[jss,iss], lon_rho[jss,iee]],
         [lat_rho[jss,iss], lat_rho[jss,iee]], '-w')
ax.plot([lon_rho[jss,iss], lon_rho[jee,iss]],
         [lat_rho[jss,iss], lat_rho[jee,iss]], '-w')
ax.plot([lon_rho[jee,iss], lon_rho[jee,iee]],
         [lat_rho[jee,iss], lat_rho[jee,iee]], '-w')
ax.plot([lon_rho[jee,iee], lon_rho[jss,iee]],
         [lat_rho[jee,iee], lat_rho[jss,iee]], '-w')


# In[20]:


fig,ax=plt.subplots(1,1)
mesh=ax.pcolormesh(lon_rho[jss:jee,iss:iee],
               lat_rho[jss:jee,iss:iee],
               h[jss:jee,iss:iee])
cbar=plt.colorbar(mesh,ax=ax)
cbar.set_label('bathy [m]')
viz_tools.plot_coastline(ax,b,coords='map')
ax.set_xlim([-125.1,-124.4])
ax.set_ylim([48.3,48.9])
ax.plot(X[:,10], Y[:,10],'g-')


# This region looks reasonable. I need to go in the mouth a little bit because our BCs cover the first 10 grid points. 
# 
# ### What's next?
# * What frequency is appropriate? We can have access to hourly snapshots.
# * Develop some scripts that would extract this region and the variables I want. Python or bash?

# ## Check out sudomain file

# In[21]:


f2 = nc.Dataset('/ocean/nsoontie/MEOPAR/LiveOcean/20160608/ocean_his_0002_UBC.nc')
f2.variables


# In[22]:


ocean_time2 = f2.variables['ocean_time']
time2 = nc.num2date(ocean_time2[:], ocean_time2.units)
time2


# In[23]:


sal2 = f2.variables['salt']
lon_rho2 = f2.variables['lon_rho']
lat_rho2 = f2.variables['lat_rho']
s_rho2 = f2.variables['s_rho']


# In[24]:


sal2


# In[25]:


sal


# In[26]:


fig,ax=plt.subplots(1,1)
mesh=ax.pcolormesh(lon_rho2, lat_rho2, sal2[0,0,:,:])
cbar=plt.colorbar(mesh,ax=ax)
cbar.set_label('salinity')
viz_tools.plot_coastline(ax,b,coords='map')
ax.set_xlim([-125,-124.4])
ax.set_ylim([48.3,48.9])
ax.plot(X[:,10], Y[:,10],'g-')


# In[27]:


f.close()


# In[28]:


UBC_subdomain.get_UBC_subdomain(['/ocean/nsoontie/MEOPAR/LiveOcean/20160608/ocean_his_0025.nc'])


# In[ ]: