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

# Experiemtne with HYCOM netcdf opening

# In[1]:


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

from salishsea_tools.nowcast import figures, residuals
from salishsea_tools import viz_tools
get_ipython().run_line_magic('matplotlib', 'inline')


# In[2]:


date = '20150730'
url = ('http://nomads.ncep.noaa.gov:9090/dods/rtofs/rtofs_global{}/rtofs_glo_2ds_forecast_3hrly_diag'.format(date))


# In[3]:


f = nc.Dataset(url)


# In[4]:


ssh=f.variables['ssh']
time = f.variables['time']
dates=nc.num2date(time[:],time.units)
lon=f.variables['lon']


# In[5]:


f


# In[6]:


lon


# In[7]:


dates


# In[8]:


print time


# In[9]:


print time[:]


# In[10]:


start = datetime.datetime(1,1,1)
dates2 = [start + datetime.timedelta(days = t) for t in time]


# In[11]:


dates2


# Ok, so there is a problem with how I was reading the date in the text file....

# In[12]:


grid_b = nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')

fig,axs= plt.subplots(1,2,figsize=(10,5))
ax=axs[0]
axm=axs[1]
date = datetime.datetime(2015,7,20)

#Hycom
iss = np.arange(1934,1920,-1)
jss = np.arange(1661,1665)
f = nc.Dataset(url)
ssh = f.variables['ssh']
time = f.variables['time']
dates = nc.num2date(time[:],time.units)
lat = f.variables['lat']
lon=f.variables['lon']
for i,j in zip(iss,jss):
    #filename = read_url(date,i,j) read website
    ax.plot(dates,ssh[:,0,j,i],label='Hycom')
    axm.plot(lon[i]-360,lat[j],'o')
viz_tools.plot_coastline(axm,grid_b,coords='map')


#Neah Bay forecast
filename = '/ocean/nsoontie/MEOPAR/sshNeahBay/txt/sshNB_2015-07-29_21.txt'
NBdata = residuals._load_surge_data(filename)
surge, dates = residuals._retrieve_surge(NBdata, datetime.datetime(2015,7,29))
ax.plot(dates[:],surge[:],label = 'Neah Bay')

# Neah Bay observations
obs = figures.get_NOAA_wlevels(figures.SITES['Neah Bay']['stn_no'], '28-Jul-2015', '31-Jul-2015')
tides = figures.get_NOAA_tides(figures.SITES['Neah Bay']['stn_no'], '28-Jul-2015', '31-Jul-2015')

res = residuals.calculate_residual(obs.wlev, obs.time, tides.pred, tides.time)
ax.plot(obs.time,res,label='observation')


axm.set_xlim([-125.5,-124])
axm.set_ylim([48,49])
ax.legend(loc=0)
ax.grid()
fig.autofmt_xdate()
ax.set_xlim([datetime.datetime(2015,7,28),datetime.datetime(2015,8,3)])


# Hycom forecast slighlty better than NOAA Jul 30/31, but everthing is small
# * July 30 forecast gives forecast for July 30 and nowards. 

# In[13]:


date = '20150731'
url = ('http://nomads.ncep.noaa.gov:9090/dods/rtofs/rtofs_global{}/rtofs_glo_2ds_forecast_3hrly_diag'.format(date))
f = nc.Dataset(url)
ssh=f.variables['ssh']
time = f.variables['time']
dates=nc.num2date(time[:],time.units)
lon=f.variables['lon']


# In[14]:


grid_b = nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')

fig,axs= plt.subplots(1,2,figsize=(10,5))
ax=axs[0]
axm=axs[1]
date = datetime.datetime(2015,7,20)

#Hycom
iss = np.arange(1934,1920,-1)
jss = np.arange(1661,1665)
f = nc.Dataset(url)
ssh = f.variables['ssh']
time = f.variables['time']
dates = nc.num2date(time[:],time.units)
lat = f.variables['lat']
lon=f.variables['lon']
for i,j in zip(iss,jss):
    #filename = read_url(date,i,j) read website
    ax.plot(dates,ssh[:,0,j,i],label='Hycom')
    axm.plot(lon[i]-360,lat[j],'o')
viz_tools.plot_coastline(axm,grid_b,coords='map')


#Neah Bay forecast
filename = '/ocean/nsoontie/MEOPAR/sshNeahBay/txt/sshNB_2015-07-30_21.txt'
NBdata = residuals._load_surge_data(filename)
surge, dates = residuals._retrieve_surge(NBdata, datetime.datetime(2015,7,29))
ax.plot(dates[:],surge[:],label = 'Neah Bay')

# Neah Bay observations
obs = figures.get_NOAA_wlevels(figures.SITES['Neah Bay']['stn_no'], '28-Jul-2015', '31-Jul-2015')
tides = figures.get_NOAA_tides(figures.SITES['Neah Bay']['stn_no'], '28-Jul-2015', '31-Jul-2015')

res = residuals.calculate_residual(obs.wlev, obs.time, tides.pred, tides.time)
ax.plot(obs.time,res,label='observation')


axm.set_xlim([-125.5,-124])
axm.set_ylim([48,49])
ax.legend(loc=0)
ax.grid()
fig.autofmt_xdate()
ax.set_xlim([datetime.datetime(2015,7,28),datetime.datetime(2015,8,3)])


# Jul 31 forecast matches well with NOAA in first part of day. Match with observations is ok...
# * July 31 forecast is for July 31 and onwards

# #HYCOM Nowcast

# In[15]:


date = '20150730'
url = ('http://nomads.ncep.noaa.gov:9090/dods/rtofs/rtofs_global{}/rtofs_glo_2ds_nowcast_3hrly_diag'.format(date))
f = nc.Dataset(url)
ssh=f.variables['ssh']
time = f.variables['time']
dates=nc.num2date(time[:],time.units)
lon=f.variables['lon']


# In[16]:


grid_b = nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')

fig,axs= plt.subplots(1,2,figsize=(10,5))
ax=axs[0]
axm=axs[1]
date = datetime.datetime(2015,7,20)

#Hycom
iss = np.arange(1934,1920,-1)
jss = np.arange(1661,1665)
f = nc.Dataset(url)
ssh = f.variables['ssh']
time = f.variables['time']
dates = nc.num2date(time[:],time.units)
lat = f.variables['lat']
lon=f.variables['lon']
for i,j in zip(iss,jss):
    #filename = read_url(date,i,j) read website
    ax.plot(dates,ssh[:,0,j,i],label='Hycom')
    axm.plot(lon[i]-360,lat[j],'o')
viz_tools.plot_coastline(axm,grid_b,coords='map')


#Neah Bay forecast
filename = '/ocean/nsoontie/MEOPAR/sshNeahBay/txt/sshNB_2015-07-30_21.txt'
NBdata = residuals._load_surge_data(filename)
surge, dates = residuals._retrieve_surge(NBdata, datetime.datetime(2015,7,29))
ax.plot(dates[:],surge[:],label = 'Neah Bay')

# Neah Bay observations
obs = figures.get_NOAA_wlevels(figures.SITES['Neah Bay']['stn_no'], '28-Jul-2015', '30-Jul-2015')
tides = figures.get_NOAA_tides(figures.SITES['Neah Bay']['stn_no'], '28-Jul-2015', '30-Jul-2015')

res = residuals.calculate_residual(obs.wlev, obs.time, tides.pred, tides.time)
ax.plot(obs.time,res,label='observation')


axm.set_xlim([-125.5,-124])
axm.set_ylim([48,49])
ax.legend(loc=0)
ax.grid()
fig.autofmt_xdate()
ax.set_xlim([datetime.datetime(2015,7,28),datetime.datetime(2015,8,3)])


# Hycom nowcast is not an observation but doesnt' look bad.
# 
# * Accessing July 30 nowcast gives data for July 28-29

# In[17]:


date = '20150731'
url = ('http://nomads.ncep.noaa.gov:9090/dods/rtofs/rtofs_global{}/rtofs_glo_2ds_nowcast_3hrly_diag'.format(date))
f = nc.Dataset(url)
ssh=f.variables['ssh']
time = f.variables['time']
dates=nc.num2date(time[:],time.units)
lon=f.variables['lon']


# In[18]:


grid_b = nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')

fig,axs= plt.subplots(1,2,figsize=(10,5))
ax=axs[0]
axm=axs[1]
date = datetime.datetime(2015,7,20)

#Hycom
iss = np.arange(1934,1920,-1)
jss = np.arange(1661,1665)
f = nc.Dataset(url)
ssh = f.variables['ssh']
time = f.variables['time']
dates = nc.num2date(time[:],time.units)
lat = f.variables['lat']
lon=f.variables['lon']
for i,j in zip(iss,jss):
    #filename = read_url(date,i,j) read website
    ax.plot(dates,ssh[:,0,j,i],label='Hycom')
    axm.plot(lon[i]-360,lat[j],'o')
viz_tools.plot_coastline(axm,grid_b,coords='map')


#Neah Bay forecast
filename = '/ocean/nsoontie/MEOPAR/sshNeahBay/txt/sshNB_2015-07-30_21.txt'
NBdata = residuals._load_surge_data(filename)
surge, dates = residuals._retrieve_surge(NBdata, datetime.datetime(2015,7,29))
ax.plot(dates[:],surge[:],label = 'Neah Bay')

# Neah Bay observations
obs = figures.get_NOAA_wlevels(figures.SITES['Neah Bay']['stn_no'], '28-Jul-2015', '30-Jul-2015')
tides = figures.get_NOAA_tides(figures.SITES['Neah Bay']['stn_no'], '28-Jul-2015', '30-Jul-2015')

res = residuals.calculate_residual(obs.wlev, obs.time, tides.pred, tides.time)
ax.plot(obs.time,res,label='observation')


axm.set_xlim([-125.5,-124])
axm.set_ylim([48,49])
ax.legend(loc=0)
ax.grid()
fig.autofmt_xdate()
ax.set_xlim([datetime.datetime(2015,7,28),datetime.datetime(2015,8,3)])


# * Accessing July 31 nowcast gives data for July 29 and 30.

# #Temperature and salinty

# In[19]:


date = '20150730'
url = ('http://nomads.ncep.noaa.gov:9090/dods/rtofs/rtofs_global{}/rtofs_glo_3dz_nowcast_6hrly_reg2'.format(date))


# In[20]:


f=nc.Dataset(url)


# In[21]:


f


# In[22]:


sal=f.variables['salinity']
lon=f.variables['lon']
lat=f.variables['lat']
dep=f.variables['lev']


# In[23]:


plt.pcolormesh(lon[:], lat[:],sal[1,0,:,:])
plt.axis([-125,-123,48,50])
plt.colorbar()


# In[24]:


iss = np.where(np.logical_and(lon[:]<-124.5, lon[:] > -125))
jss = np.where(np.logical_and(lat[:] > 48.4, lat[:] < 48.6))


# In[25]:


jss


# In[26]:


iss


# In[27]:


plt.pcolormesh(lon[:], lat[:],sal[1,0,:,:])
plt.axis([-125,-123,48,50])
plt.colorbar()

for i,j in zip(iss[0][0:-1],jss[0][::-1]+1):
    plt.plot(lon[i],lat[j],'o')


# In[28]:


fig,axs = plt.subplots(1,2,figsize=(15,5))
axs[1].pcolormesh(lon[:], lat[:],sal[1,0,:,:])
axs[1].set_xlim([-125.5,-123])
axs[1].set_ylim([48,50])

new_sal = np.ma.zeros((40,5))
new_lat=np.zeros(5)
count=0
for i,j in zip(iss[0]-1,jss[0][::-1]+1):
    new_sal[:,count] = sal[1,:,j,i]
    new_lat[count] = lat[j]
    count=count+1
    axs[1].plot(lon[i],lat[j],'o')
mesh=axs[0].pcolormesh(new_lat, dep[:],new_sal)
axs[0].set_ylim([300,0])
plt.colorbar(mesh,ax=axs[0])




# Is it possible to interpolate this product onto our bcs?

# In[29]:


for i in iss[0]:
    for j in jss[0]:
        plt.plot(sal[1,:,j,i],dep[:])
plt.axis([32,34,300,0])


# I don't think this is getting as fresh as our model is in the surface. 
# 
# Does this model account for river discharge? Probably not. 
# 
# Maybe it could be useful for the deep water an upwelling? But it isn't very deep...

# In[31]:


plt.pcolormesh(lon[:], lat[:],sal[1,0,:,:])
plt.axis([-125,-123,48,50])
plt.colorbar()

for i in iss[0]:
    for j in jss[0]:
        plt.plot(lon[i],lat[j],'o')


# #Next
# 
# * More systematic comparison of HYCOM forecasts with observations - Perhaps 3 houlry averaged obs for fairness?
# * Decide about HYCOM nowcasts - is it worth using those when we could use direct observations? But we don't have a full observation for our nowcast system.

# In[ ]: