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

# In[1]:


import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
from salishsea_tools import (tidetools, geo_tools, viz_tools)
import numpy.ma as ma
import pandas as pd
import datetime
import pytz
import os
get_ipython().run_line_magic('matplotlib', 'inline')


# In[2]:


from IPython.display import HTML

HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')


# In[2]:


grid = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/bathymetry_201702.nc')
bathy, X, Y = tidetools.get_bathy_data(grid)


# In[3]:


mesh = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/mesh_mask201702.nc')


# In[4]:


HINDCAST_PATH= '/results/SalishSea/nowcast-green/'


# In[5]:


f = pd.read_excel('/ocean/vdo/MEOPAR/2016_Nutrients_20180110_CN_edits.xlsx')
f.keys()


# In[6]:


f = f.drop(f.keys()[11:], axis=1)


# In[7]:


f.shape


# In[8]:


f = f.dropna(subset = ['Date', 'Time of Sample', 'Latitude', 'Longitude', 'Depth', 'NO3+NO'])


# In[9]:


f[:5]


# In[10]:


local = pytz.timezone ("America/Los_Angeles")


# In[11]:


import datetime


# In[12]:


datetimes = np.array([])
for index in f.index:
    dt = datetime.datetime.combine(pd.to_datetime(pd.Timestamp(f['Date'][index])), 
                          f['Time of Sample'][index])
    datetimes = np.append(datetimes, dt)


# In[13]:


f.shape


# In[14]:


datetimes.shape


# In[15]:


f = f.assign(datetime = datetimes)


# In[16]:


f.shape


# In[17]:


f.Crew.unique()


# In[63]:


list_of_lons = np.array([])
list_of_lats = np.array([])
list_of_datetimes = np.array([])
list_of_cs_ni = np.array([])
list_of_model_ni = np.array([])
list_of_depths = np.array([])
cb_lons = np.array([])
cb_lats = np.array([])
cb_datetimes = np.array([])
cb_cs_ni = np.array([])
cb_model_ni = np.array([])
cb_depths = np.array([])
for n in f.index:
    Yind, Xind = geo_tools.find_closest_model_point(f.Longitude[n], 
                                                        f.Latitude[n], 
                                                        X, Y, land_mask = bathy.mask)
    if f['Depth'][n] == 0:
        depth = 0
    elif f['Depth'][n] == 20:
        depth = 18
    if mesh.variables['tmask'][0,depth,Yind, Xind] == 1:
        date = local.localize(f['datetime'][n], is_dst=True).astimezone(pytz.utc)
        sub_dir = date.strftime('%d%b%y').lower()
        datestr = date.strftime('%Y%m%d')
        fname = 'SalishSea_1h_{}_{}_ptrc_T.nc'.format(datestr, datestr)
        nuts = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir, fname))
        if date.minute < 30:
            before = datetime.datetime(year = date.year, month = date.month, day = date.day, 
                               hour = (date.hour), minute = 30) - datetime.timedelta(hours=1)
            after = before + datetime.timedelta(hours=1)
            sub_dir2 = after.strftime('%d%b%y').lower()
            datestr2 = after.strftime('%Y%m%d')
            fname2 = 'SalishSea_1h_{}_{}_ptrc_T.nc'.format(datestr2, datestr2)
            nuts2 = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir2, fname2))
            delta = (date.minute + 30) / 60
            ni_val = ((1-delta)*(nuts.variables['nitrate'][before.hour, depth, Yind, Xind] ) + 
                       (delta)*(nuts2.variables['nitrate'][after.hour, depth, Yind, Xind] ))
        if date.minute >= 30:
            before = datetime.datetime(year = date.year, month = date.month, day = date.day, 
                                   hour = (date.hour), minute = 30)
            after = before + datetime.timedelta(hours=1)
            sub_dir2 = after.strftime('%d%b%y').lower()
            datestr2 = after.strftime('%Y%m%d')
            fname2 = 'SalishSea_1h_{}_{}_ptrc_T.nc'.format(datestr2, datestr2)
            nuts2 = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir2, fname2))
            delta = (date.minute) / 60
            ni_val = (delta*(nuts.variables['nitrate'][before.hour, depth, Yind, Xind] ) + 
                       (1- delta)*(nuts2.variables['nitrate'][after.hour, depth, Yind, Xind] ))

        list_of_lons = np.append(list_of_lons, f.Longitude[n])
        list_of_lats = np.append(list_of_lats, f.Latitude[n])
        list_of_datetimes = np.append(list_of_datetimes, date)
        if f['NO3+NO'][n] == '<0':
            list_of_cs_ni = np.append(list_of_cs_ni, 0)
        else:
            list_of_cs_ni = np.append(list_of_cs_ni, float(f['NO3+NO'][n]))
        list_of_model_ni = np.append(list_of_model_ni, ni_val)
        list_of_depths = np.append(list_of_depths, depth)
        if f.Crew[n] == 'Cowichan Bay':
            cb_lons = np.append(cb_lons, f.Longitude[n])
            cb_lats = np.append(cb_lats, f.Latitude[n])
            cb_depths = np.append(cb_depths, depth)
            cb_datetimes = np.append(cb_datetimes, date)
            if f['NO3+NO'][n] == '<0':
                cb_cs_ni = np.append(cb_cs_ni, 0)
            else:
                cb_cs_ni = np.append(cb_cs_ni, float(f['NO3+NO'][n]))
            cb_model_ni = np.append(cb_model_ni, ni_val)


# In[64]:


list_of_lats.shape


# In[65]:


list_of_model_ni.shape


# In[66]:


cb_cs_ni.shape


# In[24]:


import matplotlib as mpl
mpl.rcParams['font.size'] = 12
mpl.rcParams['axes.titlesize'] = 12


# In[95]:


fig, ax = plt.subplots(figsize = (10,10))
viz_tools.set_aspect(ax, coords = 'map')
ax.plot(list_of_lons, list_of_lats, 'ro')
viz_tools.plot_coastline(ax, grid, coords = 'map')
ax.set_ylim(48.5, 50.8)
ax.set_xlim(-125.5, -122.5);


# In[26]:


list_of_cs_ni.shape


# In[89]:


fig, ax = plt.subplots(figsize = (10,10))
ax.plot(list_of_cs_ni[list_of_depths == 18], list_of_model_ni[list_of_depths == 18],  
        'b.', alpha = 0.5, label = '20m')
ax.plot(cb_cs_ni[cb_depths == 18], cb_model_ni[cb_depths == 18],  
        '.', color = 'purple', alpha = 0.5, label = '20m CB')
ax.plot(list_of_cs_ni[list_of_depths == 0], list_of_model_ni[list_of_depths==0],  
        'g.', alpha = 0.5, label = 'surface')
ax.plot(cb_cs_ni[cb_depths == 0], cb_model_ni[cb_depths == 0],  
        '.', color = 'orange', alpha = 0.5, label = 'surface CB')
ax.plot(np.arange(0,30), color = 'grey')
ax.grid('on')
ax.set_title('Citizen Science Nitrate 2016')
ax.set_xlabel('Citizen Science')
ax.set_ylabel('Nowcast-green');
ax.legend()
print('surface bias =  ' + str(-np.mean(list_of_cs_ni[list_of_depths == 0]) 
                               + np.mean(list_of_model_ni[list_of_depths == 0])))
print('surface RMSE = ' + str(np.sqrt(np.sum((list_of_model_ni[list_of_depths == 0] 
                                              - list_of_cs_ni[list_of_depths == 0])**2) /
                              len(list_of_cs_ni[list_of_depths == 0]))))
xbar = np.mean(list_of_cs_ni[list_of_depths == 0])
print('surface Willmott = ' + str(1-(np.sum((list_of_model_ni[list_of_depths == 0] 
                                             - list_of_cs_ni[list_of_depths == 0])**2)  / 
                             np.sum((np.abs(list_of_model_ni[list_of_depths == 0] - xbar) 
                                     + np.abs(list_of_cs_ni[list_of_depths == 0] - xbar))**2))))

print('20m bias =  ' + str(-np.mean(list_of_cs_ni[list_of_depths == 18]) 
                               + np.mean(list_of_model_ni[list_of_depths == 18])))
print('20m RMSE = ' + str(np.sqrt(np.sum((list_of_model_ni[list_of_depths == 18] 
                                              - list_of_cs_ni[list_of_depths == 18])**2) /
                              len(list_of_cs_ni[list_of_depths == 18]))))
xbar = np.mean(list_of_cs_ni[list_of_depths == 18])
print('20m Willmott = ' + str(1-(np.sum((list_of_model_ni[list_of_depths == 18] 
                                             - list_of_cs_ni[list_of_depths == 18])**2)  / 
                             np.sum((np.abs(list_of_model_ni[list_of_depths == 18] - xbar) 
                                     + np.abs(list_of_cs_ni[list_of_depths == 18] - xbar))**2))))

print('bias =  ' + str(-np.mean(list_of_cs_ni) + np.mean(list_of_model_ni)))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_ni - list_of_cs_ni)**2) /
                              len(list_of_cs_ni))))
xbar = np.mean(list_of_cs_ni)
print('Willmott = ' + str(1-(np.sum((list_of_model_ni - list_of_cs_ni)**2)  / 
                             np.sum((np.abs(list_of_model_ni - xbar) 
                                     + np.abs(list_of_cs_ni - xbar))**2))))


# In[68]:


fig, ax = plt.subplots(figsize = (20,8))
ax.plot(list_of_lats[list_of_depths == 18], 
        list_of_model_ni[list_of_depths == 18] 
        - list_of_cs_ni[list_of_depths == 18], 'bo', alpha =0.5, label = '20m')
ax.plot(list_of_lats[list_of_depths == 0], 
        list_of_model_ni[list_of_depths == 0] 
        - list_of_cs_ni[list_of_depths == 0], 'go', alpha =0.5, label = 'surface')
ax.plot(cb_lats[cb_depths == 18], 
        cb_model_ni[cb_depths == 18] 
        - cb_cs_ni[cb_depths == 18], 'o', color = 'purple', alpha =0.5, label = '20m CB')
ax.plot(cb_lats[cb_depths == 0], 
        cb_model_ni[cb_depths == 0] 
        - cb_cs_ni[cb_depths == 0], 'o', color = 'orange', alpha =0.5, label = 'surface CB')
ax.legend()
ax.grid('on')
ax.set_xlabel('lat', fontsize = 15)
ax.set_ylabel('Model - Observed',fontsize = 15)
ax.set_title('Ni', fontsize = 20);


# In[69]:


fig, ax = plt.subplots(figsize = (20,8))
ax.plot(list_of_lons[list_of_depths == 18], 
        list_of_model_ni[list_of_depths == 18] 
        - list_of_cs_ni[list_of_depths == 18], 'bo', alpha =0.5, label = '20m')
ax.plot(list_of_lons[list_of_depths == 0], 
        list_of_model_ni[list_of_depths == 0] 
        - list_of_cs_ni[list_of_depths == 0], 'go', alpha =0.5, label = 'surface')
ax.plot(cb_lons[cb_depths == 18], 
        cb_model_ni[cb_depths == 18] 
        - cb_cs_ni[cb_depths == 18], 'o', color = 'purple', alpha =0.5, label = '20m CB')
ax.plot(cb_lons[cb_depths == 0], 
        cb_model_ni[cb_depths == 0] 
        - cb_cs_ni[cb_depths == 0], 'o', color = 'orange', alpha =0.5, label = 'surface CB')
ax.legend()
ax.grid('on')
ax.set_xlabel('lon', fontsize = 15)
ax.set_ylabel('Model - Observed',fontsize = 15)
ax.set_title('Ni', fontsize = 20);


# In[70]:


fig, ax = plt.subplots(figsize = (20,8))
ax.plot(list_of_datetimes[list_of_depths == 18], 
        list_of_model_ni[list_of_depths == 18] 
        - list_of_cs_ni[list_of_depths == 18], 'bo', alpha =0.5, label = '20m')
ax.plot(list_of_datetimes[list_of_depths == 0], 
        list_of_model_ni[list_of_depths == 0] 
        - list_of_cs_ni[list_of_depths == 0], 'go', alpha =0.5, label = 'surface')
ax.plot(cb_datetimes[cb_depths == 18], 
        cb_model_ni[cb_depths == 18] 
        - cb_cs_ni[cb_depths == 18], 'o', color = 'purple', alpha =0.5, label = '20m CB')
ax.plot(cb_datetimes[cb_depths == 0], 
        cb_model_ni[cb_depths == 0] 
        - cb_cs_ni[cb_depths == 0], 'o', color = 'orange', alpha =0.5, label = 'surface CB')
ax.legend()
ax.grid('on')
ax.set_xlabel('date', fontsize = 15)
ax.set_ylabel('Model - Observed',fontsize = 15)
ax.set_title('Ni', fontsize = 20);


# In[26]:


g = pd.read_excel('/ocean/vdo/MEOPAR/2016_Nutrients_20180110_CN_edits.xlsx', sheetname =1)
g.keys()


# In[27]:


g[:5]


# In[28]:


g['Site ID'].unique()


# In[29]:


g.shape


# In[30]:


g.keys()


# In[31]:


g = g.drop(g.keys()[9:], axis=1)
g.shape


# In[32]:


g = g.dropna(subset = ['DDMMYYYY', 'Time of Sample', 'Latitude', 'Longitude', 'Depth', 'SiO2 µM'])
g.shape


# In[33]:


g.Depth.unique()


# In[34]:


g.loc[g['Longitude'] < -200, 'Longitude'] = -123.31833333


# In[35]:


g.loc[g['Latitude'] < 30, 'Latitude'] = 48.7616666666667


# In[36]:


datetimes = np.array([])
for index in g.index:
    dt = datetime.datetime.combine(pd.to_datetime(pd.Timestamp(g['DDMMYYYY'][index])), 
                          g['Time of Sample'][index])
    datetimes = np.append(datetimes, dt)
g = g.assign(datetime = datetimes)


# In[37]:


g.shape


# In[71]:


list_of_lons2 = np.array([])
list_of_lats2 = np.array([])
list_of_datetimes2 = np.array([])
list_of_cs_si = np.array([])
list_of_model_si = np.array([])
list_of_depths2 = np.array([])
cb_lons2 = np.array([])
cb_lats2 = np.array([])
cb_datetimes2 = np.array([])
cb_cs_si = np.array([])
cb_model_si = np.array([])
cb_depths2 = np.array([])
for n in g.index:
    Yind, Xind = geo_tools.find_closest_model_point(g.Longitude[n], g.Latitude[n], 
                                                    X, Y, land_mask = bathy.mask)
    if g['Depth'][n] == 0:
        depth = 0
    elif g['Depth'][n] == 20:
        depth = 18
    elif g['Depth'][n] == 10:
        depth = 9
    elif g['Depth'][n] == 5:
        depth = 4
    elif g['Depth'][n] == 30:
        depth = 21
    elif g['Depth'][n] == 'S':
        depth = 0
    if mesh.variables['tmask'][0,depth,Yind, Xind] == 1:
        date = local.localize(g['datetime'][n], is_dst=True).astimezone(pytz.utc)
        sub_dir = date.strftime('%d%b%y').lower()
        datestr = date.strftime('%Y%m%d')
        fname = 'SalishSea_1h_{}_{}_ptrc_T.nc'.format(datestr, datestr)
        nuts = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir, fname))
        if date.minute < 30:
            before = datetime.datetime(year = date.year, month = date.month, day = date.day, 
                               hour = (date.hour), minute = 30) - datetime.timedelta(hours=1)
            after = before + datetime.timedelta(hours=1)
            sub_dir2 = after.strftime('%d%b%y').lower()
            datestr2 = after.strftime('%Y%m%d')
            fname2 = 'SalishSea_1h_{}_{}_ptrc_T.nc'.format(datestr2, datestr2)
            nuts2 = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir2, fname2))
            delta = (date.minute + 30) / 60
            si_val = ((1-delta)*(nuts.variables['silicon'][before.hour, depth, Yind, Xind] ) + 
                       (delta)*(nuts2.variables['silicon'][after.hour, depth, Yind, Xind] ))
        if date.minute >= 30:
            before = datetime.datetime(year = date.year, month = date.month, day = date.day, 
                                   hour = (date.hour), minute = 30)
            after = before + datetime.timedelta(hours=1)
            sub_dir2 = after.strftime('%d%b%y').lower()
            datestr2 = after.strftime('%Y%m%d')
            fname2 = 'SalishSea_1h_{}_{}_ptrc_T.nc'.format(datestr2, datestr2)
            nuts2 = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir2, fname2))
            delta = (date.minute) / 60
            si_val = (delta*(nuts.variables['silicon'][before.hour, depth, Yind, Xind] ) + 
                       (1- delta)*(nuts2.variables['silicon'][after.hour, depth, Yind, Xind] ))

        list_of_lons2 = np.append(list_of_lons2, g.Longitude[n])
        list_of_lats2 = np.append(list_of_lats2, g.Latitude[n])
        list_of_datetimes2 = np.append(list_of_datetimes2, date)
        if g['SiO2 µM'][n] < 0:
            list_of_cs_si = np.append(list_of_cs_si, 0)
        else:
            list_of_cs_si = np.append(list_of_cs_si, float(g['SiO2 µM'][n]))
        list_of_model_si = np.append(list_of_model_si, si_val)
        list_of_depths2 = np.append(list_of_depths2, depth)
        if g['Site ID'][n][:2] == 'CB':
            cb_lons2 = np.append(cb_lons2, g.Longitude[n])
            cb_lats2 = np.append(cb_lats2, g.Latitude[n])
            cb_depths2 = np.append(cb_depths2, depth)
            cb_datetimes2 = np.append(cb_datetimes2, date)
            if g['SiO2 µM'][n] <0:
                cb_cs_si = np.append(cb_cs_si, 0)
            else:
                cb_cs_si = np.append(cb_cs_si, float(g['SiO2 µM'][n]))
            cb_model_si = np.append(cb_model_si, si_val)


# In[93]:


fig, ax = plt.subplots(figsize = (10,10))
viz_tools.set_aspect(ax, coords = 'map')
ax.plot(list_of_lons2, list_of_lats2, 'ro')
viz_tools.plot_coastline(ax, grid, coords = 'map')
ax.set_title('Si stations')
ax.set_ylim(48.5, 50.8)
ax.set_xlim(-125.5, -122.5);


# In[96]:


import matplotlib.patches as patches


# In[120]:


fig, ax = plt.subplots(figsize = (10,10))
ax.plot(list_of_cs_si[list_of_depths2 == 21], list_of_model_si[list_of_depths2 == 21],  
        'o', color = 'navy', alpha = 0.5, label = '30m')
ax.plot(list_of_cs_si[list_of_depths2 == 18], list_of_model_si[list_of_depths2 == 18],  
        'b.', alpha = 0.5, label = '20m')
ax.plot(list_of_cs_si[list_of_depths2 == 9], list_of_model_si[list_of_depths2 == 9],  
        'o', color = 'yellow', alpha = 0.5, label = '10m')
ax.plot(list_of_cs_si[list_of_depths2 == 4], list_of_model_si[list_of_depths2 == 4],  
        'ro', alpha = 0.5, label = '5m')
ax.plot(list_of_cs_si[list_of_depths2 == 0], list_of_model_si[list_of_depths2==0],  
        'g.', alpha = 0.5, label = 'surface')
ax.plot(cb_cs_si[cb_depths2 == 18], cb_model_si[cb_depths2 == 18],  
        '.', color = 'purple', alpha = 0.5, label = '20m CB')
ax.plot(cb_cs_si[cb_depths2 == 0], cb_model_si[cb_depths2 == 0],  
        '.', color = 'orange', alpha = 0.5, label = 'surface CB')
ax.plot(list_of_cs_si[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19) 
                      & (list_of_depths2 == 18)], 
        list_of_model_si[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 18)],  
        '*', color = 'deeppink', alpha = 0.5, label = 'box - 20m')
ax.plot(list_of_cs_si[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19) 
                      & (list_of_depths2 == 0)], 
        list_of_model_si[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 0)],  
        '*', color = 'maroon', alpha = 0.5, label = 'box - surface')
ax.add_patch(
    patches.Rectangle(
        (22, 3),
        43,
        16,
        fill=False, linewidth=5      # remove background
    ))
ax.plot(np.arange(0,50), color = 'grey')
ax.grid('on')
ax.set_title('Citizen Science Silicon 2016')
ax.set_xlabel('Citizen Science')
ax.set_ylabel('Nowcast-green');
ax.legend()
print('surface bias =  ' + str(-np.mean(list_of_cs_si[list_of_depths2 == 0]) 
                               + np.mean(list_of_model_si[list_of_depths2 == 0])))
print('surface RMSE = ' + str(np.sqrt(np.sum((list_of_model_si[list_of_depths2 == 0] 
                                              - list_of_cs_si[list_of_depths2 == 0])**2) /
                              len(list_of_cs_si[list_of_depths2 == 0]))))
xbar = np.mean(list_of_cs_si[list_of_depths2 == 0])
print('surface Willmott = ' + str(1-(np.sum((list_of_model_si[list_of_depths2 == 0] 
                                             - list_of_cs_si[list_of_depths2 == 0])**2)  / 
                             np.sum((np.abs(list_of_model_si[list_of_depths2 == 0] - xbar) 
                                     + np.abs(list_of_cs_si[list_of_depths2 == 0] - xbar))**2))))

print('20m bias =  ' + str(-np.mean(list_of_cs_si[list_of_depths2 == 18]) 
                               + np.mean(list_of_model_si[list_of_depths2 == 18])))
print('20m RMSE = ' + str(np.sqrt(np.sum((list_of_model_si[list_of_depths2 == 18] 
                                              - list_of_cs_si[list_of_depths2 == 18])**2) /
                              len(list_of_cs_si[list_of_depths2 == 18]))))
xbar = np.mean(list_of_cs_si[list_of_depths2 == 18])
print('20m Willmott = ' + str(1-(np.sum((list_of_model_si[list_of_depths2 == 18] 
                                             - list_of_cs_si[list_of_depths2 == 18])**2)  / 
                             np.sum((np.abs(list_of_model_si[list_of_depths2 == 18] - xbar) 
                                     + np.abs(list_of_cs_si[list_of_depths2 == 18] - xbar))**2))))

print('bias =  ' + str(-np.mean(list_of_cs_si) + np.mean(list_of_model_si)))
print('RMSE = ' + str(np.sqrt(np.sum((list_of_model_si - list_of_cs_si)**2) /
                              len(list_of_cs_si))))
xbar = np.mean(list_of_cs_si)
print('Willmott = ' + str(1-(np.sum((list_of_model_si - list_of_cs_si)**2)  / 
                             np.sum((np.abs(list_of_model_si - xbar) 
                                     + np.abs(list_of_cs_si - xbar))**2))))


# In[73]:


np.histogram(list_of_depths2, bins = [-1, 3, 7, 10, 20, 25])


# # 411 samples at surface, 403 samples at 20m, 1 sample at 5, 10, and 30m

# In[121]:


fig, ax = plt.subplots(figsize = (20,8))
ax.plot(list_of_lats2[list_of_depths2 == 21], 
        list_of_model_si[list_of_depths2 == 21] 
        - list_of_cs_si[list_of_depths2 == 21], 'yo', alpha = 0.5, label = '30m')
ax.plot(list_of_lats2[list_of_depths2 == 18], 
        list_of_model_si[list_of_depths2 == 18] 
        - list_of_cs_si[list_of_depths2 == 18], 'bo', alpha =0.5, label = '20m')
ax.plot(list_of_lats2[list_of_depths2 == 9], 
        list_of_model_si[list_of_depths2 == 9] 
        - list_of_cs_si[list_of_depths2 == 9], 'o', color = 'purple', alpha =0.5, label = '10m')
ax.plot(list_of_lats2[list_of_depths2 == 4], 
        list_of_model_si[list_of_depths2 == 4] 
        - list_of_cs_si[list_of_depths2 == 4], 'ro', alpha =0.5, label = '5m')
ax.plot(list_of_lats2[list_of_depths2 == 0], 
        list_of_model_si[list_of_depths2 == 0] 
        - list_of_cs_si[list_of_depths2 == 0], 'go', alpha =0.5, label = 'surface')
ax.plot(cb_lats2[cb_depths2 == 18], 
        cb_model_si[cb_depths2 == 18] 
        - cb_cs_si[cb_depths2 == 18], 'o', color = 'purple', alpha =0.5, label = '20m CB')
ax.plot(cb_lats2[cb_depths2 == 0], 
        cb_model_si[cb_depths2 == 0] 
        - cb_cs_si[cb_depths2 == 0], 'o', color = 'orange', alpha =0.5, label = 'surface CB')
ax.plot(list_of_lats2[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 18)], 
        list_of_model_si[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 18)] 
        - list_of_cs_si[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 18)], 'o',color = 'deeppink', 
        alpha = 0.5, label = 'box - 20m')
ax.plot(list_of_lats2[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 0)], 
        list_of_model_si[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 0)] 
        - list_of_cs_si[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 0)], 'o',color = 'maroon', 
        alpha = 0.5, label = 'box - surface')
ax.legend()
ax.grid('on')
ax.set_xlabel('lat', fontsize = 15)
ax.set_ylabel('Model - Observed',fontsize = 15)
ax.set_title('Si', fontsize = 20);


# In[125]:


fig, ax = plt.subplots(figsize = (20,8))
ax.plot(list_of_lons2[list_of_depths2 == 21], 
        list_of_model_si[list_of_depths2 == 21] 
        - list_of_cs_si[list_of_depths2 == 21], 'yo', alpha = 0.5, label = '30m')
ax.plot(list_of_lons2[list_of_depths2 == 18], 
        list_of_model_si[list_of_depths2 == 18] 
        - list_of_cs_si[list_of_depths2 == 18], 'bo', alpha =0.5, label = '20m')
ax.plot(list_of_lons2[list_of_depths2 == 9], 
        list_of_model_si[list_of_depths2 == 9] 
        - list_of_cs_si[list_of_depths2 == 9], 'o', color = 'purple', alpha =0.5, label = '10m')
ax.plot(list_of_lons2[list_of_depths2 == 4], 
        list_of_model_si[list_of_depths2 == 4] 
        - list_of_cs_si[list_of_depths2 == 4], 'ro', alpha =0.5, label = '5m')
ax.plot(list_of_lons2[list_of_depths2 == 0], 
        list_of_model_si[list_of_depths2 == 0] 
        - list_of_cs_si[list_of_depths2 == 0], 'go', alpha =0.5, label = 'surface')
ax.plot(cb_lons2[cb_depths2 == 18], 
        cb_model_si[cb_depths2 == 18] 
        - cb_cs_si[cb_depths2 == 18], 'o', color = 'purple', alpha =0.5, label = '20m CB')
ax.plot(cb_lons2[cb_depths2 == 0], 
        cb_model_si[cb_depths2 == 0] 
        - cb_cs_si[cb_depths2 == 0], 'o', color = 'orange', alpha =0.5, label = 'surface CB')
ax.plot(list_of_lons2[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 18)], 
        list_of_model_si[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 18)] 
        - list_of_cs_si[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 18)], 'o',color = 'deeppink', 
        alpha = 0.5, label = 'box - 20m')
ax.plot(list_of_lons2[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 0)], 
        list_of_model_si[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 0)] 
        - list_of_cs_si[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 0)], 'o',color = 'maroon', 
        alpha = 0.5, label = 'box - surface')
ax.legend()
ax.grid('on')
ax.set_xlabel('lon', fontsize = 15)
ax.set_ylabel('Model - Observed',fontsize = 15)
ax.set_title('Si', fontsize = 20);


# In[124]:


fig, ax = plt.subplots(figsize = (20,8))
ax.plot(list_of_datetimes2[list_of_depths2 == 21], 
        list_of_model_si[list_of_depths2 == 21] 
        - list_of_cs_si[list_of_depths2 == 21], 'yo', alpha = 0.5, label = '30m')
ax.plot(list_of_datetimes2[list_of_depths2 == 18], 
        list_of_model_si[list_of_depths2 == 18] 
        - list_of_cs_si[list_of_depths2 == 18], 'bo', alpha =0.5, label = '20m')
ax.plot(list_of_datetimes2[list_of_depths2 == 9], 
        list_of_model_si[list_of_depths2 == 9] 
        - list_of_cs_si[list_of_depths2 == 9], 'o', color = 'purple', alpha =0.5, label = '10m')
ax.plot(list_of_datetimes2[list_of_depths2 == 4], 
        list_of_model_si[list_of_depths2 == 4] 
        - list_of_cs_si[list_of_depths2 == 4], 'ro', alpha =0.5, label = '5m')
ax.plot(list_of_datetimes2[list_of_depths2 == 0], 
        list_of_model_si[list_of_depths2 == 0] 
        - list_of_cs_si[list_of_depths2 == 0], 'go', alpha =0.5, label = 'surface')
ax.plot(cb_datetimes2[cb_depths2 == 18], 
        cb_model_si[cb_depths2 == 18] 
        - cb_cs_si[cb_depths2 == 18], 'o', color = 'purple', alpha =0.5, label = '20m CB')
ax.plot(cb_datetimes2[cb_depths2 == 0], 
        cb_model_si[cb_depths2 == 0] 
        - cb_cs_si[cb_depths2 == 0], 'o', color = 'orange', alpha =0.5, label = 'surface CB')
ax.plot(list_of_datetimes2[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 18)], 
        list_of_model_si[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 18)] 
        - list_of_cs_si[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 18)], 'o',color = 'deeppink', 
        alpha = 0.5, label = 'box - 20m')
ax.plot(list_of_datetimes2[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 0)], 
        list_of_model_si[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 0)] 
        - list_of_cs_si[(list_of_cs_si >= 22) & (list_of_cs_si <=65) &
                     (list_of_model_si >=3) & (list_of_model_si <= 19)
                        & (list_of_depths2 == 0)], 'o',color = 'maroon', 
        alpha = 0.5, label = 'box - surface')
ax.legend()
ax.grid('on')
ax.set_xlabel('date', fontsize = 15)
ax.set_ylabel('Model - Observed',fontsize = 15)
ax.set_title('Si', fontsize = 20);


# In[46]:


f.keys()


# In[47]:


g.keys()


# In[48]:


g.loc[g['Depth'] == 'S', 'Depth'] = 0


# In[49]:


both = pd.merge(f, g, on=['datetime', 'Latitude', 'Longitude', 'Depth']).drop_duplicates()
both.shape


# In[50]:


both[:5]


# In[77]:


list_of_cs_si3 = np.array([])
list_of_model_si3 = np.array([])
list_of_depths3 = np.array([])
list_of_cs_ni3 = np.array([])
list_of_model_ni3 = np.array([])
cb_cs_si3 = np.array([])
cb_cs_ni3 = np.array([])
cb_model_si3 = np.array([])
cb_model_ni3 = np.array([])
cb_depths3 = np.array([])
for n in both.index:
    Yind, Xind = geo_tools.find_closest_model_point(both.Longitude[n], 
                                                        both.Latitude[n], 
                                                        X, Y, land_mask = bathy.mask)
    if both['Depth'][n] == 0:
        depth = 0
    elif both['Depth'][n] == 20:
        depth = 18
    if mesh.variables['tmask'][0,depth,Yind, Xind] == 1:
        date = local.localize(both['datetime'][n], is_dst=True).astimezone(pytz.utc)
        sub_dir = date.strftime('%d%b%y').lower()
        datestr = date.strftime('%Y%m%d')
        fname = 'SalishSea_1h_{}_{}_ptrc_T.nc'.format(datestr, datestr)
        nuts = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir, fname))
        if date.minute < 30:
            before = datetime.datetime(year = date.year, month = date.month, day = date.day, 
                               hour = (date.hour), minute = 30) - datetime.timedelta(hours=1)
            after = before + datetime.timedelta(hours=1)
            sub_dir2 = after.strftime('%d%b%y').lower()
            datestr2 = after.strftime('%Y%m%d')
            fname2 = 'SalishSea_1h_{}_{}_ptrc_T.nc'.format(datestr2, datestr2)
            nuts2 = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir2, fname2))
            delta = (date.minute + 30) / 60
            si_val = ((1-delta)*(nuts.variables['silicon'][before.hour, depth, Yind, Xind] ) + 
                       (delta)*(nuts2.variables['silicon'][after.hour, depth, Yind, Xind] ))
            ni_val = ((1-delta)*(nuts.variables['nitrate'][before.hour, depth, Yind, Xind] ) + 
                       (delta)*(nuts2.variables['nitrate'][after.hour, depth, Yind, Xind] ))
        if date.minute >= 30:
            before = datetime.datetime(year = date.year, month = date.month, day = date.day, 
                                   hour = (date.hour), minute = 30)
            after = before + datetime.timedelta(hours=1)
            sub_dir2 = after.strftime('%d%b%y').lower()
            datestr2 = after.strftime('%Y%m%d')
            fname2 = 'SalishSea_1h_{}_{}_ptrc_T.nc'.format(datestr2, datestr2)
            nuts2 = nc.Dataset(os.path.join(HINDCAST_PATH, sub_dir2, fname2))
            delta = (date.minute) / 60
            si_val = (delta*(nuts.variables['silicon'][before.hour, depth, Yind, Xind] ) + 
                       (1- delta)*(nuts2.variables['silicon'][after.hour, depth, Yind, Xind] ))
            ni_val = (delta*(nuts.variables['nitrate'][before.hour, depth, Yind, Xind] ) + 
                       (1- delta)*(nuts2.variables['nitrate'][after.hour, depth, Yind, Xind] ))

        if both['SiO2 µM'][n] < 0:
            list_of_cs_si3 = np.append(list_of_cs_si3, 0)
        else:
            list_of_cs_si3 = np.append(list_of_cs_si3, float(both['SiO2 µM'][n]))
        if both['NO3+NO'][n] == '<0':
            list_of_cs_ni3 = np.append(list_of_cs_ni3, 0)
        else:
            list_of_cs_ni3 = np.append(list_of_cs_ni3, float(both['NO3+NO'][n]))
        list_of_model_si3 = np.append(list_of_model_si3, si_val)
        list_of_model_ni3 = np.append(list_of_model_ni3, ni_val)
        list_of_depths3 = np.append(list_of_depths3, depth)
        if both.Crew[n] == 'Cowichan Bay':
            if both['NO3+NO'][n] == '<0':
                cb_cs_ni3 = np.append(cb_cs_ni3, 0)
            else:
                cb_cs_ni3 = np.append(cb_cs_ni3, float(both['NO3+NO'][n]))
            if both['SiO2 µM'][n] < 0:
                cb_cs_si3 = np.append(cb_cs_si3, 0)
            else:
                cb_cs_si3 = np.append(cb_cs_si3, float(both['SiO2 µM'][n]))
            cb_model_ni3 = np.append(cb_model_ni3, ni_val)
            cb_model_si3 = np.append(cb_model_si3, si_val)
            cb_depths3 = np.append(cb_depths3, depth)


# In[78]:


list_of_cs_si3.shape


# In[79]:


list_of_depths3[list_of_depths3 == 0].shape


# In[80]:


list_of_depths3[list_of_depths3 == 18].shape


# In[81]:


list_of_depths3.shape


# In[82]:


list_of_model_ni3.shape


# In[83]:


list_of_model_si3.shape


# In[84]:


cb_cs_ni3.shape


# # 512 surface samples and 490 samples at 20m

# In[85]:


fig, ax = plt.subplots(1,2, figsize = (16,8))
ax[0].plot(list_of_cs_ni3[list_of_depths3 == 0], list_of_cs_si3[list_of_depths3 == 0], 
           'b.', alpha = 0.5)
ax[0].plot(list_of_cs_ni3[list_of_depths3 == 18], list_of_cs_si3[list_of_depths3 == 18], 
           'g.', alpha = 0.5)
ax[0].plot(cb_cs_ni3[cb_depths3 == 0], cb_cs_si3[cb_depths3 == 0], 
           '.', color = 'orange', alpha = 0.5)
ax[0].plot(cb_cs_ni3[cb_depths3 == 18], cb_cs_si3[cb_depths3 == 18], 
           '.', color = 'purple', alpha = 0.5)
ax[0].plot(np.unique(list_of_cs_ni3), 
           np.poly1d(np.polyfit(list_of_cs_ni3, list_of_cs_si3, 1))(np.unique(list_of_cs_ni3)))
x = np.arange(0,50)
#ax[0].plot(x,x, 'r-', alpha = 0.3)
ax[0].plot(x, 2*x, 'g-', alpha = 0.3)
ax[0].plot(x, 2*x+30, 'y-', alpha = 0.3)
ax[1].plot(list_of_model_ni3[list_of_depths3 == 0], 
           list_of_model_si3[list_of_depths3 == 0], 'b.', 
           alpha = 0.5, label = 'surface')
ax[1].plot(list_of_model_ni3[list_of_depths3 == 18], 
           list_of_model_si3[list_of_depths3 == 18], 'g.', 
           alpha = 0.5, label = '20m')
ax[1].plot(cb_model_ni3[cb_depths3 == 0], cb_model_si3[cb_depths3 == 0], 
           '.', color = 'orange', alpha = 0.5, label = 'surface CB')
ax[1].plot(cb_model_ni3[cb_depths3 == 18], cb_model_si3[cb_depths3 == 18], 
           '.', color = 'purple', alpha = 0.5, label = '20m CB')
ax[1].plot(np.unique(list_of_model_ni3), 
           np.poly1d(np.polyfit(list_of_model_ni3, 
                                list_of_model_si3, 1))(np.unique(list_of_model_ni3)))
x = np.arange(0,53)
#ax[1].plot(x,x, 'r-', alpha = 0.3, label = 'slope = 1')
ax[1].plot(x, 2*x, 'g-', alpha = 0.3, label = 'slope = 2')
ax[1].plot(x, 2*x+30, 'y-', alpha = 0.3, label = '2*N + 30')
ax[0].set_title('Citizen Science', fontsize = 16)
ax[1].set_title('Model', fontsize = 16)
for ax in ax:
    ax.grid('on')
    ax.set_ylabel('Si', fontsize = 14)
    ax.set_xlabel('N', fontsize = 14)
    ax.set_ylim(0,105)
    ax.set_xlim(0,40)
plt.legend();


# In[86]:


m1, b1 = np.polyfit(list_of_cs_ni3, list_of_cs_si3, 1)
print('CitSci slope = ' + str(m1))
print('CitSci y int = ' + str(b1))
m2, b2 = np.polyfit(list_of_model_ni3, list_of_model_si3, 1)
print('model slope = ' + str(m2))
print('model y int = ' + str(b2))


# In[ ]:





# In[61]:


import pickle


# In[87]:


output = open('model_si_2016.pkl', 'wb')
pickle.dump(list_of_model_si, output)
output.close()
output = open('model_ni_2016.pkl', 'wb')
pickle.dump(list_of_model_ni, output)
output.close()
output = open('cs_ni_2016.pkl', 'wb')
pickle.dump(list_of_cs_ni, output)
output.close()
output = open('cs_si_2016.pkl', 'wb')
pickle.dump(list_of_cs_si, output)
output.close()
output = open('cb_model_si_2016.pkl', 'wb')
pickle.dump(cb_model_si, output)
output.close()
output = open('cb_model_ni_2016.pkl', 'wb')
pickle.dump(cb_model_ni, output)
output.close()
output = open('cb_cs_ni_2016.pkl', 'wb')
pickle.dump(cb_cs_ni, output)
output.close()
output = open('cb_cs_si_2016.pkl', 'wb')
pickle.dump(cb_cs_si, output)
output.close()


# In[88]:


output = open('model_si_2016b.pkl', 'wb')
pickle.dump(list_of_model_si3, output)
output.close()
output = open('model_ni_2016b.pkl', 'wb')
pickle.dump(list_of_model_ni3, output)
output.close()
output = open('cs_ni_2016b.pkl', 'wb')
pickle.dump(list_of_cs_ni3, output)
output.close()
output = open('cs_si_2016b.pkl', 'wb')
pickle.dump(list_of_cs_si3, output)
output.close()
output = open('cb_model_si_2016b.pkl', 'wb')
pickle.dump(cb_model_si3, output)
output.close()
output = open('cb_model_ni_2016b.pkl', 'wb')
pickle.dump(cb_model_ni3, output)
output.close()
output = open('cb_cs_ni_2016b.pkl', 'wb')
pickle.dump(cb_cs_ni3, output)
output.close()
output = open('cb_cs_si_2016b.pkl', 'wb')
pickle.dump(cb_cs_si3, output)
output.close()


# In[ ]: