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

# In[14]:


import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import matplotlib as mpl
import netCDF4 as nc
import datetime as dt
from salishsea_tools import evaltools as et, places
import xarray as xr
import pandas as pd

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


# #### Example 1: load a time series at one location by accessing netCDF4 files stored on /results or /results2

# In[15]:


get_ipython().run_cell_magic('time', '', "start= dt.datetime(2017,3,1)\nend=dt.datetime(2017,9,30) # the code called below (evaltools.index_model_files) includes the end date \n                          # in the values returned\nbasedir='/results/SalishSea/nowcast-green.201812/'\nnam_fmt='nowcast'\nflen=1 # files contain 1 day of data each\nftype= 'ptrc_T' # load bio files\ntres=24 # 1: hourly resolution; 24: daily resolution   <- try changing to 1 and loading hourly data\nflist=et.index_model_files(start,end,basedir,nam_fmt,flen,ftype,tres)\n# flist contains paths: file pathes; t_0 timestemp of start of each file; t_n: timestamp of start of next file\nprint(flist)\n")


# In[16]:


# get model i,j of location S3 from places
ij,ii=places.PLACES['S3']['NEMO grid ji']
ik=0 # choose surface level


# In[17]:


bio=xr.open_mfdataset(flist['paths'])


# In[18]:


bio


# In[19]:


get_ipython().run_cell_magic('time', '', 'tt=bio.time_centered\nmicZ=bio.microzooplankton.isel(deptht=ik,y=ij,x=ii)\ndiat=bio.diatoms.isel(deptht=ik,y=ij,x=ii)\n')


# In[20]:


get_ipython().run_cell_magic('time', '', "fig,ax=plt.subplots(1,1,figsize=(12,3))\nax.plot(tt,diat,'c-',label='diatoms')\nax.plot(tt,micZ,'-',color='darkorange',label='microzooplankton')\nax.legend(loc=2);\nax.set_ylabel('Concentration ($\\mu$M N)')\nax.set_xlim(tt[0],tt[-1])\n")


# In[21]:


# Adjust date format display
fig,ax=plt.subplots(1,1,figsize=(12,3))
ax.plot(tt,diat,'c-',label='diatoms')
ax.plot(tt,micZ,'-',color='darkorange',label='microzooplankton')
ax.legend(loc=2);
ax.set_ylabel('Concentration ($\mu$M N)')
ax.set_xlim(tt[0],tt[-1])
yearsFmt = mdates.DateFormatter('%d %b')
ax.xaxis.set_major_formatter(yearsFmt)
ax.set_title('March 2017')


# In[22]:


bio.close()


# ### repeat with hourly data over a shorter interval

# In[11]:


start= dt.datetime(2017,3,1)
end=dt.datetime(2017,3,5) # the code called below (evaltools.index_model_files) includes the end date 
                          # in the values returned
tres=1 # 1: hourly resolution; 24: daily resolution   <- try changing to 1 and loading hourly data
flist=et.index_model_files(start,end,basedir,nam_fmt,flen,ftype,tres)
# flist contains paths: file pathes; t_0 timestemp of start of each file; t_n: timestamp of start of next file
flist['paths'][0]


# In[12]:


get_ipython().run_cell_magic('time', '', "bio=xr.open_mfdataset(flist['paths'])\ntt=bio.time_centered\nmicZ=bio.microzooplankton.isel(deptht=0,y=ij,x=ii)\ndiat=bio.diatoms.isel(deptht=0,y=ij,x=ii)\n")


# In[13]:


get_ipython().run_cell_magic('time', '', "fig,ax=plt.subplots(1,1,figsize=(12,3))\nax.plot(tt,diat,'c-',label='diatoms')\nax.plot(tt,micZ,'-',color='darkorange',label='microzooplankton')\nax.legend(loc=2);\nax.set_ylabel('Concentration ($\\mu$M N)')\nax.set_xlim(tt[0],tt[-1])\n")


# In[14]:


bio.close()


# ### Aside: this can also be done using ERDDAP, but only to access 201905 model output

# In[15]:


get_ipython().run_cell_magic('time', '', "hc1905 = xr.open_dataset('https://salishsea.eos.ubc.ca/erddap/griddap/ubcSSg3DBiologyFields1hV19-05')\n")


# In[16]:


get_ipython().run_cell_magic('time', '', 'hc1905\n')


# In[17]:


get_ipython().run_cell_magic('time', '', 'modt=hc1905.time.sel(time=slice(start,end))\nmicZ1905=hc1905.microzooplankton.sel(time=slice(start,end)).isel(depth=ik,gridY=ij,gridX=ii)\ndiat1905=hc1905.diatoms.sel(time=slice(start,end)).isel(depth=ik,gridY=ij,gridX=ii)\n')


# In[18]:


get_ipython().run_cell_magic('time', '', "fig,ax=plt.subplots(1,1,figsize=(12,3))\nax.plot(modt,diat1905,'c-',label='diatoms')\nax.plot(modt,micZ1905,'-',color='darkorange',label='microzooplankton')\nax.legend(loc=2);\nax.set_ylabel('Concentration ($\\mu$M N)')\nax.set_xlim(modt[0],modt[-1]);\n")


# In[19]:


hc1905.close()


# ## back to accessing files from /results and /results2, but now switch to loading 201905 results

# In[20]:


basedir='/results2/SalishSea/nowcast-green.201905/'


# ## Plot mid Feb -- Jun Surface phytoplankton and nitrate at S3 for 2017:

# In[21]:


newstart=dt.datetime(2017,2,15)
newend=dt.datetime(2017,6,1)
flist=et.index_model_files(newstart,newend,basedir,'nowcast',1,'ptrc_T',1)
bio2=xr.open_mfdataset(flist['paths'])
ij,ii=places.PLACES['S3']['NEMO grid ji']
ik=0 # choose surface level


# In[22]:


get_ipython().run_cell_magic('time', '', "fig,ax=plt.subplots(1,1,figsize=(12,3))\np1=ax.plot(bio2.time_centered,bio2.diatoms.isel(deptht=ik,y=ij,x=ii)+\\\n        bio2.flagellates.isel(deptht=ik,y=ij,x=ii)+bio2.ciliates.isel(deptht=ik,y=ij,x=ii),\n           '-',color='teal',label='Phytoplankton')\np2=ax.plot(bio2.time_centered,bio2.nitrate.isel(deptht=ik,y=ij,x=ii),\n           '-',color='orange',label='Nitrate')\nax.legend(handles=[p1[0],p2[0]],loc=1)\nax.set_ylabel('$\\mu$M N')\n")


# ## Fraser River flow for same times

# In[23]:


dfFra=pd.read_csv('/ocean/eolson/MEOPAR/obs/ECRivers/Flow/FraserHopeDaily__Dec-2-2020_10_31_05PM.csv',
                  skiprows=1)
# the original file contains both flow and water level information in the same field (Value)
# keep only the flow data, where PARAM=1 (drop PARAM=2 values)
# flow units are m3/s
dfFra.drop(dfFra.loc[dfFra.PARAM==2].index,inplace=True)  


# In[25]:


# rename 'Value' column to 'Flow' now that we have removed all the water level rows
dfFra.rename(columns={'Value':'Flow'},inplace=True)


# In[26]:


dfFra


# In[27]:


# no time information so use dt.date
dfFra['Date']=[dt.date(iyr,1,1)+dt.timedelta(days=idd-1) for iyr, idd in zip(dfFra['YEAR'],dfFra['DD'])]


# In[28]:


dfFra.head(2)


# In[29]:


# select portion of dataframe in desired date range
dfFra2=dfFra.loc[(dfFra.Date>=newstart.date())&(dfFra.Date<=newend.date())]


# In[33]:


fig,ax=plt.subplots(1,1,figsize=(12,3))
ax.plot(dfFra2['Date'],dfFra2['Flow'],'c-')
ax.set_ylabel('Flow (m$^3$s$^{-1}$)')
ax.set_title('Fraser Flow at Hope')


# In[ ]: