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

# In[1]:


import numpy as np
import matplotlib.pyplot as plt
from salishsea_tools import viz_tools, places, visualisations
import netCDF4 as nc # unless you prefer xarray
import datetime as dt
import os
import glob
import cmocean

from IPython.display import Markdown, display
get_ipython().run_line_magic('matplotlib', 'inline')

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[ ]:





# ### Load December files from the 201905 hindcast

# In[2]:


d0=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20071201_20071231.nc') #Use December 2007 for both 2007 and 2008 "winter values 


# In[3]:


d1=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20081201_20081231.nc')


# In[4]:


d2=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20091201_20091231.nc')


# In[5]:


d3=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20101201_20101231.nc')


# In[6]:


d4=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20111201_20111231.nc')


# In[7]:


d5=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20121201_20121231.nc')


# In[8]:


d6=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20131201_20131231.nc')


# In[9]:


d7=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20141201_20141231.nc')


# In[10]:


d8=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20151201_20151231.nc')


# In[11]:


d9=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20161201_20161231.nc')


# In[12]:


d10=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20171201_20171231.nc')


# In[13]:


d11=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20181201_20181231.nc')


# In[14]:


d12=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20191201_20191231.nc')


# In[15]:


d13=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20201201_20201231.nc')


# In[16]:


print(d1.variables.keys())


# In[17]:


AllYears_N_Dec=((d0.variables['nitrate'][0,...])+\
            (d0.variables['nitrate'][0,...])+\
            (d1.variables['nitrate'][0,...])+\
            (d2.variables['nitrate'][0,...])+\
            (d3.variables['nitrate'][0,...])+\
            (d4.variables['nitrate'][0,...])+\
            (d5.variables['nitrate'][0,...])+\
            (d6.variables['nitrate'][0,...])+\
            (d7.variables['nitrate'][0,...])+\
            (d8.variables['nitrate'][0,...])+\
            (d9.variables['nitrate'][0,...])+\
            (d10.variables['nitrate'][0,...])+\
            (d11.variables['nitrate'][0,...])+\
            (d12.variables['nitrate'][0,...]))/14

### use 2007 values for both 2007 and 2008 because there are no December 2006 files


# In[18]:


ColdYears_N_Dec=((d0.variables['nitrate'][0,...])+\
              (d2.variables['nitrate'][0,...])+\
              (d3.variables['nitrate'][0,...])+\
              (d4.variables['nitrate'][0,...]))/4


# In[19]:


WarmYears_N_Dec=((d8.variables['nitrate'][0,...])+\
              (d11.variables['nitrate'][0,...])+\
              (d12.variables['nitrate'][0,...])+\
              (d13.variables['nitrate'][0,...]))/4


# In[20]:


AllYears_Si_Dec=((d0.variables['silicon'][0,...])+\
                 (d0.variables['silicon'][0,...])+\
            (d1.variables['silicon'][0,...])+\
            (d2.variables['silicon'][0,...])+\
            (d3.variables['silicon'][0,...])+\
            (d4.variables['silicon'][0,...])+\
            (d5.variables['silicon'][0,...])+\
            (d6.variables['silicon'][0,...])+\
            (d7.variables['silicon'][0,...])+\
            (d8.variables['silicon'][0,...])+\
            (d9.variables['silicon'][0,...])+\
            (d10.variables['silicon'][0,...])+\
            (d11.variables['silicon'][0,...])+\
            (d12.variables['silicon'][0,...]))/14


# In[21]:


ColdYears_Si_Dec=((d0.variables['silicon'][0,...])+\
              (d2.variables['silicon'][0,...])+\
              (d3.variables['silicon'][0,...])+\
              (d4.variables['silicon'][0,...]))/4


# In[22]:


WarmYears_Si_Dec=((d8.variables['silicon'][0,...])+\
              (d11.variables['silicon'][0,...])+\
              (d12.variables['silicon'][0,...])+\
              (d13.variables['silicon'][0,...]))/4


# ### Load January files from the 201905 hindcast

# In[23]:


j0=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20070101_20070131.nc')


# In[24]:


j1=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20080101_20080131.nc')


# In[25]:


j2=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20090101_20090131.nc')


# In[26]:


j3=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20100101_20100131.nc')


# In[27]:


j4=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20110101_20110131.nc')


# In[28]:


j5=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20120101_20120131.nc')


# In[29]:


j6=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20130101_20130131.nc')


# In[30]:


j7=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20140101_20140131.nc')


# In[31]:


j8=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20150101_20150131.nc')


# In[32]:


j9=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20160101_20160131.nc')


# In[33]:


j10=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20170101_20170131.nc')


# In[34]:


j11=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20180101_20180131.nc')


# In[35]:


j12=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20190101_20190131.nc')


# In[36]:


j13=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20200101_20200131.nc')


# In[37]:


AllYears_N_Jan=((j0.variables['nitrate'][0,...])+\
            (j1.variables['nitrate'][0,...])+\
            (j2.variables['nitrate'][0,...])+\
            (j3.variables['nitrate'][0,...])+\
            (j4.variables['nitrate'][0,...])+\
            (j5.variables['nitrate'][0,...])+\
            (j6.variables['nitrate'][0,...])+\
            (j7.variables['nitrate'][0,...])+\
            (j8.variables['nitrate'][0,...])+\
            (j9.variables['nitrate'][0,...])+\
            (j10.variables['nitrate'][0,...])+\
            (j11.variables['nitrate'][0,...])+\
            (j12.variables['nitrate'][0,...])+\
            (j13.variables['nitrate'][0,...]))/14


# In[38]:


ColdYears_N_Jan=((j1.variables['nitrate'][0,...])+\
              (j3.variables['nitrate'][0,...])+\
              (j4.variables['nitrate'][0,...])+\
              (j5.variables['nitrate'][0,...]))/4


# In[39]:


WarmYears_N_Jan=((j8.variables['nitrate'][0,...])+\
              (j11.variables['nitrate'][0,...])+\
              (j12.variables['nitrate'][0,...])+\
              (j13.variables['nitrate'][0,...]))/4


# In[40]:


AllYears_Si_Jan=((j0.variables['silicon'][0,...])+\
            (j1.variables['silicon'][0,...])+\
            (j2.variables['silicon'][0,...])+\
            (j3.variables['silicon'][0,...])+\
            (j4.variables['silicon'][0,...])+\
            (j5.variables['silicon'][0,...])+\
            (j6.variables['silicon'][0,...])+\
            (j7.variables['silicon'][0,...])+\
            (j8.variables['silicon'][0,...])+\
            (j9.variables['silicon'][0,...])+\
            (j10.variables['silicon'][0,...])+\
            (j11.variables['silicon'][0,...])+\
            (j12.variables['silicon'][0,...])+\
            (j13.variables['silicon'][0,...]))/14


# In[41]:


ColdYears_Si_Jan=((j1.variables['silicon'][0,...])+\
              (j3.variables['silicon'][0,...])+\
              (j4.variables['silicon'][0,...])+\
              (j5.variables['silicon'][0,...]))/4


# In[42]:


WarmYears_Si_Jan=((j8.variables['silicon'][0,...])+\
              (j11.variables['silicon'][0,...])+\
              (j12.variables['silicon'][0,...])+\
              (j13.variables['silicon'][0,...]))/4


# ### Load February files from the 201905 hindcast

# In[43]:


f0=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20070201_20070228.nc')


# In[44]:


f1=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20080201_20080229.nc')


# In[45]:


f2=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20090201_20090228.nc')


# In[46]:


f3=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20100201_20100228.nc')


# In[47]:


f4=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20110201_20110228.nc')


# In[48]:


f5=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20120201_20120229.nc')


# In[49]:


f6=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20130201_20130228.nc')


# In[50]:


f7=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20140201_20140228.nc')


# In[51]:


f8=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20150201_20150228.nc')


# In[52]:


f9=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20160201_20160229.nc')


# In[53]:


f10=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20170201_20170228.nc')


# In[54]:


f11=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20180201_20180228.nc')


# In[55]:


f12=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20190201_20190228.nc')


# In[56]:


f13=nc.Dataset('/results2/SalishSea/month-avg.201905/SalishSeaCast_1m_ptrc_T_20200201_20200229.nc')


# In[57]:


print(f1.variables.keys())


# In[58]:


AllYears_N_Feb=((f0.variables['nitrate'][0,...])+\
            (f1.variables['nitrate'][0,...])+\
            (f2.variables['nitrate'][0,...])+\
            (f3.variables['nitrate'][0,...])+\
            (f4.variables['nitrate'][0,...])+\
            (f5.variables['nitrate'][0,...])+\
            (f6.variables['nitrate'][0,...])+\
            (f7.variables['nitrate'][0,...])+\
            (f8.variables['nitrate'][0,...])+\
            (f9.variables['nitrate'][0,...])+\
            (f10.variables['nitrate'][0,...])+\
            (f11.variables['nitrate'][0,...])+\
            (f12.variables['nitrate'][0,...])+\
             (f13.variables['nitrate'][0,...])   )/14

### use 2007 values for both 2007 and 2008 because there are no December 2006 files


# In[59]:


ColdYears_N_Feb=((f1.variables['nitrate'][0,...])+\
              (f3.variables['nitrate'][0,...])+\
              (f4.variables['nitrate'][0,...])+\
              (f5.variables['nitrate'][0,...]))/4


# In[60]:


WarmYears_N_Feb=((f8.variables['nitrate'][0,...])+\
              (f11.variables['nitrate'][0,...])+\
              (f12.variables['nitrate'][0,...])+\
              (f13.variables['nitrate'][0,...]))/4


# In[61]:


AllYears_Si_Feb=((f0.variables['silicon'][0,...])+\
            (f1.variables['silicon'][0,...])+\
            (f2.variables['silicon'][0,...])+\
            (f3.variables['silicon'][0,...])+\
            (f4.variables['silicon'][0,...])+\
            (f5.variables['silicon'][0,...])+\
            (f6.variables['silicon'][0,...])+\
            (f7.variables['silicon'][0,...])+\
            (f8.variables['silicon'][0,...])+\
            (f9.variables['silicon'][0,...])+\
            (f10.variables['silicon'][0,...])+\
            (f11.variables['silicon'][0,...])+\
            (f12.variables['silicon'][0,...])+\
            (f13.variables['silicon'][0,...]))/14


# In[62]:


ColdYears_Si_Feb=((f1.variables['silicon'][0,...])+\
              (f3.variables['silicon'][0,...])+\
              (f4.variables['silicon'][0,...])+\
              (f5.variables['silicon'][0,...]))/4


# In[63]:


WarmYears_Si_Feb=((f8.variables['silicon'][0,...])+\
              (f11.variables['silicon'][0,...])+\
              (f12.variables['silicon'][0,...])+\
              (f13.variables['silicon'][0,...]))/4


# In[64]:


AllYears_N_Winter=(AllYears_N_Dec+AllYears_N_Jan+AllYears_N_Feb)/3


# In[65]:


AllYears_Si_Winter=(AllYears_Si_Dec+AllYears_Si_Jan+AllYears_Si_Feb)/3


# In[66]:


ColdYears_N_Winter=(ColdYears_N_Dec+ColdYears_N_Jan+ColdYears_N_Feb)/3


# In[67]:


WarmYears_N_Winter=(WarmYears_N_Dec+WarmYears_N_Jan+WarmYears_N_Feb)/3


# In[68]:


ColdYears_Si_Winter=(ColdYears_Si_Dec+ColdYears_Si_Jan+ColdYears_Si_Feb)/3


# In[69]:


WarmYears_Si_Winter=(WarmYears_Si_Dec+WarmYears_Si_Jan+WarmYears_Si_Feb)/3


# ### Thalweg plot
#   method using contour_thalweg from visualisations.py in tools repo

# In[70]:


#open bathy file and meshmask
fbathy=nc.Dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/bathymetry_201702.nc')
fmesh=nc.Dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc')


# In[71]:


#il=0
fig,ax=plt.subplots(1,1,figsize=(11,3))
#fig.title('Mesozooplankton Biomass')
levels = np.arange(-2, 2, .01)


cb=visualisations.contour_thalweg(ax,((ColdYears_N_Winter)-(AllYears_N_Winter)),fbathy,fmesh,clevels=levels,cmap=cmocean.cm.balance)
cb.set_label('$\mu$M', rotation=90,labelpad=6)

fig.suptitle('Winter Nitrate Anomalies - Cold Years', fontsize=16)
#ax.text(15, 100, '(d)', fontsize=15,  color='k')
#plt.savefig('MesozoopThalweg_Fig6d.png', bbox_inches='tight',dpi=1000,transparent=False)


# In[72]:


#il=0
fig,ax=plt.subplots(1,1,figsize=(11,3))
#fig.title('Mesozooplankton Biomass')
levels = np.arange(-2.5, 2.5, .01)


cb=visualisations.contour_thalweg(ax,((WarmYears_N_Winter)-(AllYears_N_Winter)),fbathy,fmesh,clevels=levels,cmap=cmocean.cm.balance)
cb.set_label('$\mu$M', rotation=90,labelpad=6)

fig.suptitle('Winter Nitrate Anomalies - Warm Years', fontsize=16)
#ax.text(15, 100, '(d)', fontsize=15,  color='k')
#plt.savefig('MesozoopThalweg_Fig6d.png', bbox_inches='tight',dpi=1000,transparent=False)


# In[73]:


#il=0
fig,ax=plt.subplots(1,1,figsize=(11,3))
#fig.title('Mesozooplankton Biomass')
levels = np.arange(-3.5, 3.5, .05)


cb=visualisations.contour_thalweg(ax,((ColdYears_Si_Winter)-(AllYears_Si_Winter)),fbathy,fmesh,clevels=levels,cmap=cmocean.cm.balance)
cb.set_label('$\mu$M', rotation=90,labelpad=6)

fig.suptitle('Winter Silicon Anomalies - Cold Years', fontsize=16)
#ax.text(15, 100, '(d)', fontsize=15,  color='k')
#plt.savefig('MesozoopThalweg_Fig6d.png', bbox_inches='tight',dpi=1000,transparent=False)


# In[74]:


#il=0
fig,ax=plt.subplots(1,1,figsize=(11,3))
#fig.title('Mesozooplankton Biomass')
levels = np.arange(-3.5, 3.5, .05)


cb=visualisations.contour_thalweg(ax,((WarmYears_Si_Winter)-(AllYears_Si_Winter)),fbathy,fmesh,clevels=levels,cmap=cmocean.cm.balance)
cb.set_label('$\mu$M', rotation=90,labelpad=6)

fig.suptitle('Winter Silicon Anomalies - Warm Years', fontsize=16)
#ax.text(15, 100, '(d)', fontsize=15,  color='k')
#plt.savefig('MesozoopThalweg_Fig6d.png', bbox_inches='tight',dpi=1000,transparent=False)


# In[98]:


fig,ax=plt.subplots(4,1,figsize=(11,15))
#fig.title('Mesozooplankton Biomass')
levels = np.arange(-2.5, 2.5, .01)

fig.tight_layout(pad=2.5)

cb=visualisations.contour_thalweg(ax[0],((ColdYears_N_Winter)-(AllYears_N_Winter)),fbathy,fmesh,clevels=levels,cmap=cmocean.cm.balance)
cb.set_label('Nitrate - Cold Years ($\mu$M)', rotation=90,labelpad=6)
#fig.suptitle('Winter Nitrate Anomalies - Cold Years', fontsize=16)
ax[0].text(15, -10, '(a)', fontsize=15,  color='k')


#fig,ax=plt.subplots(1,1,figsize=(11,3))
levels = np.arange(-2.5, 2.5, .01)
cb=visualisations.contour_thalweg(ax[1],((WarmYears_N_Winter)-(AllYears_N_Winter)),fbathy,fmesh,clevels=levels,cmap=cmocean.cm.balance)
cb.set_label('Nitrate - Warm Years ($\mu$M)', rotation=90,labelpad=6)
#fig.suptitle('Winter Nitrate Anomalies - Warm Years', fontsize=16)
ax[1].text(15, -10, '(b)', fontsize=15,  color='k')


#fig,ax=plt.subplots(1,1,figsize=(11,3))
levels = np.arange(-3.5, 3.5, .05)
cb=visualisations.contour_thalweg(ax[2],((ColdYears_Si_Winter)-(AllYears_Si_Winter)),fbathy,fmesh,clevels=levels,cmap=cmocean.cm.balance)
cb.set_label('Silicon - Cold Years ($\mu$M)', rotation=90,labelpad=6)
#fig.suptitle('Winter Silicon Anomalies - Cold Years', fontsize=16)
ax[2].text(15, -10, '(c)', fontsize=15,  color='k')


#fig,ax=plt.subplots(1,1,figsize=(11,3))
levels = np.arange(-3.5, 3.5, .05)
cb=visualisations.contour_thalweg(ax[3],((WarmYears_Si_Winter)-(AllYears_Si_Winter)),fbathy,fmesh,clevels=levels,cmap=cmocean.cm.balance)
cb.set_label('Silicon - Warm Years ($\mu$M)', rotation=90,labelpad=6)
#ax[3].suptitle('Winter Silicon Anomalies - Warm Years', fontsize=16)
ax[3].text(15, -10, '(d)', fontsize=15,  color='k')

#plt.savefig('SuppFigureS5_NutrientThalwegPlots.png', bbox_inches='tight',dpi=1000,transparent=False)


# In[86]:


f0.close()


# In[87]:


f1.close()


# In[88]:


f2.close()


# In[89]:


f3.close()


# In[90]:


f4.close()


# In[91]:


f5.close()


# In[92]:


f6.close()


# In[93]:


f7.close()


# In[94]:


f8.close()


# In[95]:


f9.close()


# In[96]:


f10.close()


# In[97]:


f11.close()


# In[98]:


f12.close()


# In[99]:


f13.close()


# In[ ]: