import xarray as xr
import datetime as dt
import cmocean
import dateutil as dutil
from salishsea_tools import viz_tools, places, evaltools as et
import matplotlib.dates as mdates
import numpy as np
from matplotlib import pyplot as plt
import matplotlib as mpl
import datetime as dt
mpl.rc('xtick', labelsize=14)
mpl.rc('ytick', labelsize=14)
mpl.rc('legend', fontsize=16)
mpl.rc('axes', titlesize=16)
mpl.rc('figure', titlesize=16)
mpl.rc('axes', labelsize=16)
mpl.rc('font', size=16)
mpl.rc('legend', numpoints= 1)
mpl.rc('lines', markersize= 8)
%matplotlib inline
def newlight(mumax,a,Imax,I):
#a=alpha/mumax
alpha=a*mumax
beta=alpha/(np.exp(Imax*alpha/mumax)-1)
b=beta/mumax
f=(alpha+beta)/alpha*((alpha+beta)/beta)**(beta/alpha)
mu=mumax*(1-np.exp(-a*I))*(np.exp(-b*I))*f
return mu
I=np.arange(0,700)
mumaxH={'d':5.486816753e-05,'f':1.945325940e-05,'d2':6.0349996E-05}
IoptH={'d':62.658287255,'f':14.918639823,'d2':67.13387920}
aH={'d':0.07215007215,'f':0.30303030303,'d2':0.067340067}
mumaxJ={'d':5.7e-05,'f':1.8e-05}
IoptJ={'d':75,'f':20}
aJ={'d':0.052,'f':0.25}
fig,ax=plt.subplots(1,1,figsize=(16,6))
s='d'
ax.plot(I,newlight(mumaxH[s],aH[s],IoptH[s],I)*24*3600,'b-',label='diatoms H')
ax.axvline(IoptH[s],color='r',linestyle=':')
#s='d2'
#ax.plot(I,newlight(mumaxH[s],aH[s],IoptH[s],I)*24*3600,'b--',label='diatoms 1905')
#ax.axvline(IoptH[s],color='r',linestyle=':')
s='f'
ax.plot(I,newlight(mumaxH[s],aH[s],IoptH[s],I)*24*3600,'r-',label='flags H')
ax.axvline(IoptH[s],color='b',linestyle=':')
s='d'
ax.plot(I,newlight(mumaxJ[s],aJ[s],IoptJ[s],I)*24*3600,'c-',label='diatoms new')
ax.axvline(IoptJ[s],color='c',linestyle=':')
#s='f'
#ax.plot(I,newlight(mumaxJ[s],aJ[s],IoptJ[s],I)*24*3600,'m-',label='flags J')
#ax.axvline(IoptJ[s],color='m',linestyle=':')
ax.legend()
ax.set_xlabel('W/m$^2$')
ax.set_ylabel('1/d')
ax.set_title('Light-dependent growth')
Text(0.5, 1.0, 'Light-dependent growth')
fig,ax=plt.subplots(1,1,figsize=(16,6))
s='d'
ax.plot(I,newlight(mumaxH[s],aH[s],IoptH[s],I)*24*3600,'b-',label='diatoms H')
ax.axvline(IoptH[s],color='r',linestyle=':')
s='f'
ax.plot(I,newlight(mumaxH[s],aH[s],IoptH[s],I)*24*3600,'r-',label='flags H')
ax.axvline(IoptH[s],color='b',linestyle=':')
s='d'
ax.plot(I,newlight(mumaxJ[s],aJ[s],IoptJ[s],I)*24*3600,'c-',label='diatoms J')
ax.axvline(IoptJ[s],color='c',linestyle=':')
s='f'
ax.plot(I,newlight(mumaxJ[s],aJ[s],IoptJ[s],I)*24*3600,'m-',label='flags J')
ax.axvline(IoptJ[s],color='m',linestyle=':')
ax.legend()
ax.set_xlim(0,80)
ax.set_xlabel('W/m$^2$')
ax.set_ylabel('1/d')
ax.set_title('Light-dependent growth')
Text(0.5, 1.0, 'Light-dependent growth')
ij,ii=places.PLACES['S3']['NEMO grid ji']
rundir='/results2/SalishSea/nowcast-green.201905/'
startdate=dt.datetime(2015,1,1)
enddate=dt.datetime(2015,7,1)
fin=et.index_model_files(startdate,enddate,rundir,'nowcast',1,'carp_T',1)
fobj=xr.open_mfdataset(fin['paths'])
tt=fobj.time_centered if 'time_centered' in fobj.coords else fobj.time_counter
par=fobj.PAR.isel(deptht=slice(0,3),y=ij,x=ii).compute()
fobj.close()
fig,ax=plt.subplots(1,1,figsize=(18,2))
ax.plot(tt,par,'k-')
[<matplotlib.lines.Line2D at 0x7f2a321d7b80>, <matplotlib.lines.Line2D at 0x7f2a3247da90>, <matplotlib.lines.Line2D at 0x7f2a3247d7c0>]
np.shape(lightold),np.shape(lightoldDay)
((960,), (40,))
len(lightoldDay),len(lightold),int(len(lightold/24))
(40, 960, 960)
mumaxJ={'d':5.8e-05,'f':1.8e-05}
IoptJ={'d':75,'f':20}
aJ={'d':0.045,'f':0.25}
mumaxJ={'d':5.8e-05,'f':1.8e-05}
IoptJ={'d':75,'f':20}
aJ={'d':0.045,'f':0.25}
fig,ax=plt.subplots(1,1,figsize=(18,3))
s='d'
ax.plot(I,newlight(mumaxH[s],aH[s],IoptH[s],I)*24*3600,'b-',label='diatoms H')
ax.axvline(IoptH[s],color='b',linestyle=':')
ax.plot(I,newlight(5.65e-05,0.052,75,I)*24*3600,'k-',label='diatoms J')
ax.axvline(75,color='k',linestyle=':')
s='d'
ax.plot(I,newlight(mumaxJ[s],aJ[s],IoptJ[s],I)*24*3600,'c-',label='diatoms J')
ax.axvline(IoptJ[s],color='c',linestyle=':')
ax.legend()
ax.set_xlim(0,600)
ax.set_xlabel('W/m$^2$')
ax.set_ylabel('1/d')
ax.set_title('Light-dependent growth')
Text(0.5, 1.0, 'Light-dependent growth')
s='d'
lightold=np.mean(newlight(mumaxH[s],aH[s],IoptH[s],par)*24*3600,1)
lightnew=np.mean(newlight(mumaxJ[s],aJ[s],IoptJ[s],par)*24*3600,1)
lightoldDay=np.array([np.mean(lightold[(5+i0*24):(5+(i0+1)*24)]) for i0 in range(0,int(len(lightold)/24))])
lightnewDay=np.array([np.mean(lightnew[(5+i0*24):(5+(i0+1)*24)]) for i0 in range(0,int(len(lightnew)/24))])
fig,ax=plt.subplots(3,1,figsize=(16,9))
s='d'
ax[0].plot(tt,lightold,'b-')
ax[0].plot(tt,lightnew,'c-')
ax[1].plot(tt[12::24],lightoldDay,'b-')
ax[1].plot(tt[12::24],lightnewDay,'c-')
ax[2].plot(tt,par,'k-')
for iax in ax:
iax.axvline(dt.datetime(2015,3,7),color='k')
iax.axvline(dt.datetime(2015,3,10),color='k')
iax.axvline(dt.datetime(2015,3,18),color='k')