import numpy as np
import matplotlib.pyplot as plt
import os
import pandas as pd
import netCDF4 as nc
import datetime as dt
from salishsea_tools import evaltools as et, viz_tools
import gsw
import matplotlib.gridspec as gridspec
import matplotlib as mpl
import matplotlib.dates as mdates
import cmocean as cmo
import scipy.interpolate as sinterp
import pickle
import cmocean
import json
import f90nml
from collections import OrderedDict
fs=16
mpl.rc('xtick', labelsize=fs)
mpl.rc('ytick', labelsize=fs)
mpl.rc('legend', fontsize=fs)
mpl.rc('axes', titlesize=fs)
mpl.rc('axes', labelsize=fs)
mpl.rc('figure', titlesize=fs)
mpl.rc('font', size=fs)
mpl.rc('text', usetex=True)
mpl.rc('text.latex', preamble = r'''
\usepackage{txfonts}
\usepackage{lmodern}
''')
mpl.rc('font', family='sans-serif', weight='normal', style='normal')
import warnings
warnings.filterwarnings('ignore')
from IPython.display import Markdown, display
%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>''')
PATH= '/results2/SalishSea/nowcast-green.201905/'
year=2007
# Parameters
year = 2017
display(Markdown('''# Year: '''+ str(year)))
start_date = dt.datetime(year,1,1)
end_date = dt.datetime(year,12,31)
flen=1
namfmt='nowcast'
filemap={'nitrate':'ptrc_T','silicon':'ptrc_T','ammonium':'ptrc_T','diatoms':'ptrc_T',
'ciliates':'ptrc_T','flagellates':'ptrc_T','vosaline':'grid_T','votemper':'grid_T'}
fdict={'ptrc_T':1,'grid_T':1}
df1=et.loadDFO(datelims=(start_date,end_date))
print(len(df1),'data points')
df1[['Year','Month','Day','Lat','Lon','Pressure','Depth','N','Si','Chlorophyll_Extracted',
'ConsT','AbsSal']].head()
1884 data points
Year | Month | Day | Lat | Lon | Pressure | Depth | N | Si | Chlorophyll_Extracted | ConsT | AbsSal | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2017.0 | 2.0 | 19.0 | 48.5 | -124.735667 | 1.0 | 1.0 | 9.52 | 26.96 | 1.12 | 8.521459 | 29.283488 |
1 | 2017.0 | 2.0 | 19.0 | 48.5 | -124.735667 | 5.6 | 5.6 | 9.53 | 26.88 | NaN | 8.525728 | 29.298387 |
2 | 2017.0 | 2.0 | 19.0 | 48.5 | -124.735667 | 9.9 | 9.9 | 9.70 | 26.49 | 0.95 | 8.535950 | 29.333996 |
3 | 2017.0 | 2.0 | 19.0 | 48.5 | -124.735667 | 20.1 | 19.9 | 11.34 | 25.35 | 0.64 | 8.475745 | 30.465823 |
4 | 2017.0 | 2.0 | 19.0 | 48.5 | -124.735667 | 30.0 | 29.7 | 12.35 | 26.16 | NaN | 8.401070 | 30.628528 |
data=et.matchData(df1,filemap,fdict,start_date,end_date,'nowcast',PATH,1,quiet=True);
# density calculations:
data['rho']=gsw.rho(data['AbsSal'],data['ConsT'],data['Pressure'])
data['mod_rho']=gsw.rho(data['mod_vosaline'],data['mod_votemper'],
gsw.p_from_z(-1*data['Z'],data['Lat']))
# load chl to N ratio from namelist
nml=f90nml.read(os.path.join(PATH,'01jan'+str(year)[-2:],'namelist_smelt_cfg'))
mod_chl_N=nml['nampisopt']['zzn2chl']
print('Parameter values from 01jan'+str(year)[-2:]+' namelist_smelt_cfg:')
print(' Chl:N = ',mod_chl_N)
print(' zz_bfsi = ',nml['nampisrem']['zz_bfsi'])
print(' zz_remin_d_bsi = ',nml['nampisrem']['zz_remin_d_bsi'])
print(' zz_w_sink_d_bsi = ',nml['nampissink']['zz_w_sink_d_bsi'])
print(' zz_alpha_b_si = ',nml['nampissink']['zz_alpha_b_si'])
print(' zz_alpha_b_d = ',nml['nampissink']['zz_alpha_b_d'])
Parameter values from 01jan17 namelist_smelt_cfg: Chl:N = 2.0 zz_bfsi = 6e-05 zz_remin_d_bsi = 1.1e-06 zz_w_sink_d_bsi = 0.00028 zz_alpha_b_si = 0.92 zz_alpha_b_d = 0.0
# chlorophyll calculations
data['l10_obsChl']=np.log10(data['Chlorophyll_Extracted']+0.01)
data['l10_modChl']=np.log10(mod_chl_N*(data['mod_diatoms']+data['mod_ciliates']+data['mod_flagellates'])+0.01)
data['mod_Chl']=mod_chl_N*(data['mod_diatoms']+data['mod_ciliates']+data['mod_flagellates'])
data['Chl']=data['Chlorophyll_Extracted']
# prep and load dictionary to save stats in
if os.path.isfile('vET-HC1905-DFO-NutChlPhys-stats.json'):
with open('vET-HC1905-DFO-NutChlPhys-stats.json', 'r') as fstat:
statsDict = json.load(fstat);
statsDict[year]=dict();
else:
statsDict={year:dict()};
cm1=cmocean.cm.thermal
theta=-30
lon0=-123.9
lat0=49.3
with nc.Dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/bathymetry_201702.nc') as bathy:
bathylon=np.copy(bathy.variables['nav_lon'][:,:])
bathylat=np.copy(bathy.variables['nav_lat'][:,:])
bathyZ=np.copy(bathy.variables['Bathymetry'][:,:])
def byDepth(ax,obsvar,modvar,lims):
ps=et.varvarPlot(ax,data,obsvar,modvar,'Z',(15,22),'z','m',('mediumseagreen','darkturquoise','navy'))
l=ax.legend(handles=ps)
ax.set_xlabel('Obs')
ax.set_ylabel('Model')
ax.plot(lims,lims,'k-',alpha=.5)
ax.set_xlim(lims)
ax.set_ylim(lims)
ax.set_aspect(1)
return ps,l
def byRegion(ax,obsvar,modvar,lims):
ps1=et.varvarPlot(ax,dJDF,obsvar,modvar,cols=('b',),lname='SJDF')
ps2=et.varvarPlot(ax,dSJGI,obsvar,modvar,cols=('c',),lname='SJGI')
ps3=et.varvarPlot(ax,dSOG,obsvar,modvar,cols=('y',),lname='SOG')
ps4=et.varvarPlot(ax,dNSOG,obsvar,modvar,cols=('m',),lname='NSOG')
l=ax.legend(handles=[ps1[0][0],ps2[0][0],ps3[0][0],ps4[0][0]])
ax.set_xlabel('Obs')
ax.set_ylabel('Model')
ax.plot(lims,lims,'k-',alpha=.5)
ax.set_xlim(lims)
ax.set_ylim(lims)
ax.set_aspect(1)
return (ps1,ps2,ps3,ps4),l
def bySeason(ax,obsvar,modvar,lims):
for axi in ax:
axi.plot(lims,lims,'k-')
axi.set_xlim(lims)
axi.set_ylim(lims)
axi.set_aspect(1)
axi.set_xlabel('Obs')
axi.set_ylabel('Model')
ps=et.varvarPlot(ax[0],JFM,obsvar,modvar,cols=('crimson','darkturquoise','navy'))
ax[0].set_title('Jan-Mar')
ps=et.varvarPlot(ax[1],Apr,obsvar,modvar,cols=('crimson','darkturquoise','navy'))
ax[1].set_title('Apr')
ps=et.varvarPlot(ax[2],MJJA,obsvar,modvar,cols=('crimson','darkturquoise','navy'))
ax[2].set_title('May-Aug')
ps=et.varvarPlot(ax[3],SOND,obsvar,modvar,cols=('crimson','darkturquoise','navy'))
ax[3].set_title('Sep-Dec')
return
def ErrErr(fig,ax,obsvar1,modvar1,obsvar2,modvar2,lims1,lims2):
m=ax.scatter(data[modvar1]-data[obsvar1],data[modvar2]-data[obsvar2],c=data['Z'],s=1,cmap='gnuplot')
cb=fig.colorbar(m,ax=ax,label='Depth (m)')
ax.set_xlim(lims1)
ax.set_ylim(lims2)
ax.set_aspect((lims1[1]-lims1[0])/(lims2[1]-lims2[0]))
return m,cb
fig, ax = plt.subplots(1,2,figsize = (13,6))
viz_tools.set_aspect(ax[0], coords = 'map')
ax[0].plot(data['Lon'], data['Lat'], 'ro',label='data')
ax[0].plot(data.loc[data.Si>75,['Lon']],data.loc[data.Si>75,['Lat']],'*',color='y',label='high Si')
grid = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/bathymetry_201702.nc')
viz_tools.plot_coastline(ax[0], grid, coords = 'map',isobath=.1)
ax[0].set_ylim(48, 50.5)
ax[0].legend()
ax[0].set_xlim(-125.7, -122.5);
ax[0].set_title('Observation Locations');
viz_tools.set_aspect(ax[1], coords = 'map')
#ax[1].plot(data['Lon'], data['Lat'], 'ro',label='data')
dJDF=data.loc[(data.Lon<-123.6)&(data.Lat<48.6)]
ax[1].plot(dJDF['Lon'],dJDF['Lat'],'b.',label='JDF')
dSJGI=data.loc[(data.Lon>=-123.6)&(data.Lat<48.9)]
ax[1].plot(dSJGI['Lon'],dSJGI['Lat'],'c.',label='SJGI')
dSOG=data.loc[(data.Lat>=48.9)&(data.Lon>-124.0)]
ax[1].plot(dSOG['Lon'],dSOG['Lat'],'y.',label='SOG')
dNSOG=data.loc[(data.Lat>=48.9)&(data.Lon<=-124.0)]
ax[1].plot(dNSOG['Lon'],dNSOG['Lat'],'m.',label='NSOG')
grid = nc.Dataset('/data/vdo/MEOPAR/NEMO-forcing/grid/bathymetry_201702.nc')
viz_tools.plot_coastline(ax[1], grid, coords = 'map')
ax[1].set_ylim(48, 50.5)
ax[1].legend()
ax[1].set_xlim(-125.7, -122.5);
# Also set up seasonal groupings:
iz=(data.Z<15)
JFM=data.loc[iz&(data.dtUTC<=dt.datetime(year,4,1)),:]
Apr=data.loc[iz&(data.dtUTC<=dt.datetime(year,5,1))&(data.dtUTC>dt.datetime(year,4,1)),:]
MJJA=data.loc[iz&(data.dtUTC<=dt.datetime(year,9,1))&(data.dtUTC>dt.datetime(year,5,1)),:]
SOND=data.loc[iz&(data.dtUTC>dt.datetime(year,9,1)),:]
statsubs=OrderedDict({'z < 15 m':data.loc[data.Z<15],
'15 m < z < 22 m':data.loc[(data.Z>=15)&(data.Z<22)],
'z >= 22 m':data.loc[data.Z>=22],
'z > 50 m':data.loc[data.Z>50],
'all':data,
'z < 15 m, JFM':JFM,
'z < 15 m, Apr':Apr,
'z < 15 m, MJJA':MJJA,
'z < 15 m, SOND': SOND})
obsvar='N'
modvar='mod_nitrate'
statsDict[year]['NO3']=OrderedDict()
for isub in statsubs:
statsDict[year]['NO3'][isub]=dict()
var=statsDict[year]['NO3'][isub]
var['N'],mmean,omean,var['Bias'],var['RMSE'],var['WSS']=et.stats(statsubs[isub].loc[:,[obsvar]],
statsubs[isub].loc[:,[modvar]])
tbl,tdf=et.displayStats(statsDict[year]['NO3'],level='Subset',suborder=list(statsubs.keys()))
tbl
Bias | N | RMSE | WSS | ||
---|---|---|---|---|---|
Subset | |||||
0 | z < 15 m | 1.21132 | 508 | 6.09019 | 0.837433 |
1 | 15 m < z < 22 m | -0.211189 | 137 | 3.84811 | 0.840336 |
2 | z >= 22 m | -0.739887 | 945 | 2.60958 | 0.893272 |
3 | z > 50 m | -1.26612 | 579 | 2.20042 | 0.892088 |
4 | all | -0.0709275 | 1590 | 4.1441 | 0.923395 |
5 | z < 15 m, JFM | -2.65867 | 42 | 3.66444 | 0.822522 |
6 | z < 15 m, Apr | 5.87382 | 185 | 7.35724 | 0.642506 |
7 | z < 15 m, MJJA | -1.51093 | 184 | 4.70592 | 0.878824 |
8 | z < 15 m, SOND | -0.841547 | 97 | 6.57179 | 0.786077 |
fig, ax = plt.subplots(1,2,figsize = (16,7))
ps,l=byDepth(ax[0],obsvar,modvar,(0,40))
ax[0].set_title('NO$_3$ ($\mu$M) By Depth')
ps,l=byRegion(ax[1],obsvar,modvar,(0,40))
ax[1].set_title('NO$_3$ ($\mu$M) By Region');
fig, ax = plt.subplots(1,4,figsize = (16,3.3))
bySeason(ax,obsvar,modvar,(0,30))
fig,ax=plt.subplots(1,1,figsize=(20,.3))
ax.plot(data.dtUTC,np.ones(np.shape(data.dtUTC)),'k.')
ax.set_xlim((dt.datetime(year,1,1),dt.datetime(year,12,31)))
ax.set_title('Data Timing')
ax.yaxis.set_visible(False)
fig,ax=plt.subplots(1,2,figsize=(12,4))
ax[0].set_xlabel('Density Error (kg m$^{-3}$)')
ax[0].set_ylabel('NO$_3$ ($\mu$M) Error')
m,cb=ErrErr(fig,ax[0],'rho','mod_rho',obsvar,modvar,(-3,3),(-15,15))
ax[1].set_xlabel('Salinity Error (g kg$^{-1}$)')
ax[1].set_ylabel('NO$_3$ ($\mu$M) Error')
m,cb=ErrErr(fig,ax[1],'AbsSal','mod_vosaline',obsvar,modvar,(-2.5,2.5),(-15,15))
obsvar='Si'
modvar='mod_silicon'
statsDict[year]['dSi']=OrderedDict()
for isub in statsubs:
statsDict[year]['dSi'][isub]=dict()
var=statsDict[year]['dSi'][isub]
var['N'],mmean,omean,var['Bias'],var['RMSE'],var['WSS']=et.stats(statsubs[isub].loc[:,[obsvar]],
statsubs[isub].loc[:,[modvar]])
tbl,tdf=et.displayStats(statsDict[year]['dSi'],level='Subset',suborder=list(statsubs.keys()))
tbl
Bias | N | RMSE | WSS | ||
---|---|---|---|---|---|
Subset | |||||
0 | z < 15 m | 1.42277 | 508 | 13.4813 | 0.6459 |
1 | 15 m < z < 22 m | -3.897 | 137 | 10.3804 | 0.684622 |
2 | z >= 22 m | -5.30834 | 944 | 8.63213 | 0.776742 |
3 | z > 50 m | -5.8247 | 578 | 8.18832 | 0.823395 |
4 | all | -3.03473 | 1589 | 10.567 | 0.781749 |
5 | z < 15 m, JFM | -8.49306 | 42 | 9.46056 | 0.771748 |
6 | z < 15 m, Apr | 14.7136 | 185 | 18.1407 | 0.478543 |
7 | z < 15 m, MJJA | -6.8016 | 184 | 10.3377 | 0.691526 |
8 | z < 15 m, SOND | -4.0315 | 97 | 9.09493 | 0.502836 |
mv=(0,80)
fig, ax = plt.subplots(1,2,figsize = (16,7))
ps,l=byDepth(ax[0],obsvar,modvar,mv)
ax[0].set_title('Dissolved Silica ($\mu$M) By Depth')
ps,l=byRegion(ax[1],obsvar,modvar,mv)
ax[1].set_title('Dissolved Silica ($\mu$M) By Region');
fig, ax = plt.subplots(1,4,figsize = (16,3.3))
bySeason(ax,obsvar,modvar,mv)
fig,ax=plt.subplots(1,1,figsize=(20,.3))
ax.plot(data.dtUTC,np.ones(np.shape(data.dtUTC)),'k.')
ax.set_xlim((dt.datetime(year,1,1),dt.datetime(year,12,31)))
ax.set_title('Data Timing')
ax.yaxis.set_visible(False)
fig,ax=plt.subplots(1,2,figsize=(12,4))
ax[0].set_xlabel('Density Error (kg m$^{-3}$)')
ax[0].set_ylabel('dSi Error ($\mu$M)')
m,cb=ErrErr(fig,ax[0],'rho','mod_rho',obsvar,modvar,(-3,3),(-25,25))
ax[1].set_xlabel('Salinity Error (g kg$^{-1}$)')
ax[1].set_ylabel('dSi Error ($\mu$M)')
m,cb=ErrErr(fig,ax[1],'AbsSal','mod_vosaline',obsvar,modvar,(-2.5,2.5),(-25,25))
fig, ax = plt.subplots(1,2,figsize = (15,8))
cols=('crimson','red','orangered','darkorange','gold','chartreuse','green','lightseagreen','cyan',
'darkturquoise','royalblue','lightskyblue','blue','darkblue','mediumslateblue','blueviolet',
'darkmagenta','fuchsia','deeppink','pink')
ii0=start_date
for ii in range(0,int((end_date-start_date).days/30)):
iii=(data.dtUTC>=(start_date+dt.timedelta(days=ii*30)))&(data.dtUTC<(start_date+dt.timedelta(days=(ii+1)*30)))
ax[0].plot(data.loc[iii,['mod_nitrate']].values-data.loc[iii,['N']].values, data.loc[iii,['Z']].values,
'.', color = cols[ii],label=str(ii))
ax[1].plot(data.loc[iii,['mod_silicon']].values-data.loc[iii,['Si']].values, data.loc[iii,['Z']].values,
'.', color = cols[ii],label=str(ii))
for axi in (ax[0],ax[1]):
axi.legend(loc=4)
axi.set_ylim(400,0)
axi.set_ylabel('Depth (m)')
ax[0].set_xlabel('Model - Obs')
ax[1].set_xlabel('Model - Obs')
ax[0].set_xlim(-15,15)
ax[1].set_xlim(-40,20)
ax[0].set_title('NO3')
ax[1].set_title('dSi')
Text(0.5, 1.0, 'dSi')
fig,ax=plt.subplots(1,2,figsize=(15,6))
p1=ax[0].plot(dJDF['N'],dJDF['Si'],'b.',label='SJDF')
p2=ax[0].plot(dSJGI['N'],dSJGI['Si'],'c.',label='SJGI')
p3=ax[0].plot(dSOG['N'],dSOG['Si'],'y.',label='SOG')
p4=ax[0].plot(dNSOG['N'],dNSOG['Si'],'m.',label='NSOG')
ax[0].plot(np.arange(0,41),1.35*np.arange(0,41)+6.46,'k-',label='OBC')
ax[0].set_title('Observed')
ax[0].set_xlabel('NO3')
ax[0].set_ylabel('dSi')
ax[0].set_xlim(0,40)
ax[0].set_ylim(0,85)
ax[0].legend()
p5=ax[1].plot(dJDF['mod_nitrate'],dJDF['mod_silicon'],'b.',label='SJDF')
p6=ax[1].plot(dSJGI['mod_nitrate'],dSJGI['mod_silicon'],'c.',label='SJGI')
p7=ax[1].plot(dSOG['mod_nitrate'],dSOG['mod_silicon'],'y.',label='SOG')
p8=ax[1].plot(dNSOG['mod_nitrate'],dNSOG['mod_silicon'],'m.',label='NSOG')
ax[1].plot(np.arange(0,41),1.35*np.arange(0,41)+6.46,'k-',label='OBC')
ax[1].set_title('Model')
ax[1].set_xlabel('NO3')
ax[1].set_ylabel('dSi')
ax[1].set_xlim(0,40)
ax[1].set_ylim(0,85)
ax[1].legend()
#ax[0].plot(np.arange(0,35),1.3*np.arange(0,35),'k-')
#ax[1].plot(np.arange(0,35),1.3*np.arange(0,35),'k-')
<matplotlib.legend.Legend at 0x7fb2d58c3520>
fig,ax=plt.subplots(1,2,figsize=(15,6))
p1=ax[0].plot(dJDF['AbsSal'], dJDF['Si']-1.3*dJDF['N'],'b.',label='SJDF')
p2=ax[0].plot(dSJGI['AbsSal'],dSJGI['Si']-1.3*dSJGI['N'],'c.',label='SJGI')
p3=ax[0].plot(dSOG['AbsSal'],dSOG['Si']-1.3*dSOG['N'],'y.',label='SOG')
p4=ax[0].plot(dNSOG['AbsSal'],dNSOG['Si']-1.3*dNSOG['N'],'m.',label='NSOG')
ax[0].set_title('Observed')
ax[0].set_xlabel('S (g/kg)')
ax[0].set_ylabel('dSi-1.3NO3')
ax[0].set_xlim(10,35)
ax[0].set_ylim(0,45)
ax[0].legend()
p5=ax[1].plot(dJDF['mod_vosaline'],dJDF['mod_silicon']-1.3*dJDF['mod_nitrate'],'b.',label='SJDF')
p6=ax[1].plot(dSJGI['mod_vosaline'],dSJGI['mod_silicon']-1.3*dSJGI['mod_nitrate'],'c.',label='SJGI')
p7=ax[1].plot(dSOG['mod_vosaline'],dSOG['mod_silicon']-1.3*dSOG['mod_nitrate'],'y.',label='SOG')
p8=ax[1].plot(dNSOG['mod_vosaline'],dNSOG['mod_silicon']-1.3*dNSOG['mod_nitrate'],'m.',label='NSOG')
ax[1].set_title('Model')
ax[1].set_xlabel('S (g/kg)')
ax[1].set_ylabel('dSi-1.3NO3')
ax[1].set_xlim(10,35)
ax[1].set_ylim(0,45)
ax[1].legend()
<matplotlib.legend.Legend at 0x7fb2ba3e7e20>
obsvar='l10_obsChl'
modvar='l10_modChl'
statsDict[year]['Chl log10']=OrderedDict()
for isub in statsubs:
statsDict[year]['Chl log10'][isub]=dict()
var=statsDict[year]['Chl log10'][isub]
var['N'],mmean,omean,var['Bias'],var['RMSE'],var['WSS']=et.stats(statsubs[isub].loc[:,[obsvar]],
statsubs[isub].loc[:,[modvar]])
obsvar='Chlorophyll_Extracted'
modvar='mod_Chl'
statsDict[year]['Chl']=OrderedDict()
for isub in statsubs:
statsDict[year]['Chl'][isub]=dict()
var=statsDict[year]['Chl'][isub]
var['N'],mmean,omean,var['Bias'],var['RMSE'],var['WSS']=et.stats(statsubs[isub].loc[:,[obsvar]],
statsubs[isub].loc[:,[modvar]])
tempD={'Chl log10':statsDict[year]['Chl log10'],'Chl':statsDict[year]['Chl']}
tbl,tdf=et.displayStatsFlex(tempD,('Variable','Subset','Metric',''),
['Order','Subset','Metric'],
['Variable','Metric'],
suborder=list(statsubs.keys()))
tbl
Variable | Chl | Chl log10 | |||||||
---|---|---|---|---|---|---|---|---|---|
Bias | N | RMSE | WSS | Bias | N | RMSE | WSS | ||
Subset | |||||||||
0 | z < 15 m | -1.54052 | 309 | 5.62303 | 0.41881 | -0.116013 | 309 | 0.566884 | 0.514615 |
1 | 15 m < z < 22 m | -0.518307 | 116 | 1.78706 | 0.424334 | -0.0796708 | 116 | 0.479648 | 0.543346 |
2 | z >= 22 m | -1.88707 | 3 | 2.57095 | 0.534625 | -0.441195 | 3 | 0.492631 | 0.611506 |
3 | z > 50 m | nan | 0 | nan | nan | nan | 0 | nan | nan |
4 | all | -1.2659 | 428 | 4.87229 | 0.452938 | -0.108443 | 428 | 0.544116 | 0.617788 |
5 | z < 15 m, JFM | -2.79809 | 40 | 7.79091 | 0.29813 | -0.309659 | 40 | 0.625575 | 0.386615 |
6 | z < 15 m, Apr | -5.43267 | 64 | 8.35052 | 0.553279 | -0.514023 | 64 | 0.695956 | 0.6122 |
7 | z < 15 m, MJJA | 0.521975 | 126 | 4.23393 | 0.36672 | 0.200752 | 126 | 0.475027 | 0.411375 |
8 | z < 15 m, SOND | -1.04019 | 79 | 2.80295 | 0.403642 | -0.200747 | 79 | 0.553642 | 0.408978 |
fig, ax = plt.subplots(1,2,figsize = (14,6))
ax[0].plot(np.arange(-.6,1.6,.1),np.arange(-.6,1.6,.1),'k-')
ps=et.varvarPlot(ax[0],data,'l10_obsChl','l10_modChl','Z',(5,10,15,20,25),'z','m',('crimson','darkorange','lime','mediumseagreen','darkturquoise','navy'))
ax[0].legend(handles=ps)
ax[0].set_xlabel('Obs')
ax[0].set_ylabel('Model')
ax[0].set_title('log10[Chl ($\mu$g/L)+0.01] By Depth')
ax[1].plot(np.arange(0,35),np.arange(0,35),'k-')
ps=et.varvarPlot(ax[1],data,'Chlorophyll_Extracted','mod_Chl','Z',(5,10,15,20,25),'z','m',('crimson','darkorange','lime','mediumseagreen','darkturquoise','navy'))
ax[1].legend(handles=ps)
ax[1].set_xlabel('Obs')
ax[1].set_ylabel('Model')
ax[1].set_title('Chl ($\mu$g/L) By Depth');
fig, ax = plt.subplots(1,2,figsize = (14,6))
obsvar='l10_obsChl'; modvar='l10_modChl'
ps,l=byRegion(ax[0],obsvar,modvar,(-.6,1.6))
ax[0].set_title('Log10 Chl ($\mu$g/L) By Region');
obsvar='Chlorophyll_Extracted'; modvar='mod_Chl'
ps,l=byRegion(ax[1],obsvar,modvar,(0,30))
ax[1].set_title('Chl ($\mu$g/L) By Region');
obsvar='ConsT'
modvar='mod_votemper'
statsDict[year]['Temperature']=OrderedDict()
for isub in statsubs:
statsDict[year]['Temperature'][isub]=dict()
var=statsDict[year]['Temperature'][isub]
var['N'],mmean,omean,var['Bias'],var['RMSE'],var['WSS']=et.stats(statsubs[isub].loc[:,[obsvar]],
statsubs[isub].loc[:,[modvar]])
tbl,tdf=et.displayStats(statsDict[year]['Temperature'],level='Subset',suborder=list(statsubs.keys()))
tbl
Bias | N | RMSE | WSS | ||
---|---|---|---|---|---|
Subset | |||||
0 | z < 15 m | 0.153269 | 655 | 1.10507 | 0.968584 |
1 | 15 m < z < 22 m | -0.0160375 | 158 | 0.511579 | 0.961782 |
2 | z >= 22 m | -0.104113 | 1022 | 0.368898 | 0.96031 |
3 | z > 50 m | -0.078967 | 647 | 0.354609 | 0.954128 |
4 | all | -0.00465754 | 1835 | 0.730905 | 0.975418 |
5 | z < 15 m, JFM | -0.520092 | 54 | 0.700092 | 0.549084 |
6 | z < 15 m, Apr | -0.00468493 | 219 | 0.416628 | 0.842328 |
7 | z < 15 m, MJJA | 0.429239 | 274 | 1.50409 | 0.911308 |
8 | z < 15 m, SOND | 0.110097 | 108 | 1.03424 | 0.932098 |
fig, ax = plt.subplots(1,2,figsize = (16,7))
ps,l=byDepth(ax[0],obsvar,modvar,(5,20))
ax[0].set_title('$\Theta$ ($^{\circ}$C) By Depth')
ps,l=byRegion(ax[1],obsvar,modvar,(5,20))
ax[1].set_title('$\Theta$ ($^{\circ}$C) By Region');
fig, ax = plt.subplots(1,4,figsize = (16,3.3))
bySeason(ax,obsvar,modvar,mv)
fig,ax=plt.subplots(1,1,figsize=(20,.3))
ax.plot(data.dtUTC,np.ones(np.shape(data.dtUTC)),'k.')
ax.set_xlim((dt.datetime(year,1,1),dt.datetime(year,12,31)))
ax.set_title('Data Timing')
ax.yaxis.set_visible(False)
obsvar='AbsSal'
modvar='mod_vosaline'
statsDict[year]['Salinity']=OrderedDict()
for isub in statsubs:
statsDict[year]['Salinity'][isub]=dict()
var=statsDict[year]['Salinity'][isub]
var['N'],mmean,omean,var['Bias'],var['RMSE'],var['WSS']=et.stats(statsubs[isub].loc[:,[obsvar]],
statsubs[isub].loc[:,[modvar]])
tbl,tdf=et.displayStats(statsDict[year]['Salinity'],level='Subset',suborder=list(statsubs.keys()))
tbl
Bias | N | RMSE | WSS | ||
---|---|---|---|---|---|
Subset | |||||
0 | z < 15 m | -0.0904083 | 655 | 1.5314 | 0.931232 |
1 | 15 m < z < 22 m | 0.0580848 | 158 | 0.530179 | 0.94627 |
2 | z >= 22 m | 0.132994 | 1022 | 0.28146 | 0.990731 |
3 | z > 50 m | 0.147509 | 647 | 0.258548 | 0.989591 |
4 | all | 0.0468009 | 1835 | 0.951543 | 0.967239 |
5 | z < 15 m, JFM | -0.50658 | 54 | 1.58329 | 0.827336 |
6 | z < 15 m, Apr | -0.13077 | 219 | 0.881454 | 0.935913 |
7 | z < 15 m, MJJA | 0.0577092 | 274 | 2.00616 | 0.924753 |
8 | z < 15 m, SOND | -0.176257 | 108 | 1.08788 | 0.891446 |
fig, ax = plt.subplots(1,2,figsize = (16,7))
ps,l=byDepth(ax[0],obsvar,modvar,(0,36))
ax[0].set_title('S$_A$ (g kg$^{-1}$) By Depth')
ps,l=byRegion(ax[1],obsvar,modvar,(0,36))
ax[1].set_title('S$_A$ (g kg$^{-1}$) By Region');
fig, ax = plt.subplots(1,4,figsize = (16,3.3))
bySeason(ax,obsvar,modvar,(0,36))
fig,ax=plt.subplots(1,1,figsize=(20,.3))
ax.plot(data.dtUTC,np.ones(np.shape(data.dtUTC)),'k.')
ax.set_xlim((dt.datetime(year,1,1),dt.datetime(year,12,31)))
ax.set_title('Data Timing')
ax.yaxis.set_visible(False)
obsvar='rho'
modvar='mod_rho'
statsDict[year]['Density']=OrderedDict()
for isub in statsubs:
statsDict[year]['Density'][isub]=dict()
var=statsDict[year]['Density'][isub]
var['N'],mmean,omean,var['Bias'],var['RMSE'],var['WSS']=et.stats(statsubs[isub].loc[:,[obsvar]],
statsubs[isub].loc[:,[modvar]])
tbl,tdf=et.displayStats(statsDict[year]['Density'],level='Subset',suborder=list(statsubs.keys()))
tbl
Bias | N | RMSE | WSS | ||
---|---|---|---|---|---|
Subset | |||||
0 | z < 15 m | -0.10353 | 655 | 1.27226 | 0.938817 |
1 | 15 m < z < 22 m | 0.0443576 | 158 | 0.454994 | 0.940907 |
2 | z >= 22 m | 0.118125 | 1022 | 0.234376 | 0.993265 |
3 | z > 50 m | 0.12618 | 647 | 0.209889 | 0.991867 |
4 | all | 0.032654 | 1835 | 0.791324 | 0.97542 |
5 | z < 15 m, JFM | -0.333153 | 54 | 1.18491 | 0.835796 |
6 | z < 15 m, Apr | -0.100339 | 219 | 0.721694 | 0.933834 |
7 | z < 15 m, MJJA | -0.0390939 | 274 | 1.68556 | 0.926471 |
8 | z < 15 m, SOND | -0.158668 | 108 | 0.922291 | 0.917526 |
fig, ax = plt.subplots(1,2,figsize = (16,7))
ps,l=byDepth(ax[0],obsvar,modvar,(1010,1030))
ax[0].set_title('Density (kg m$^{-3}$) By Depth')
ps,l=byRegion(ax[1],obsvar,modvar,(1010,1030))
ax[1].set_title('Density (kg m$^{-3}$) By Region');
fig, ax = plt.subplots(1,4,figsize = (16,3.3))
bySeason(ax,obsvar,modvar,(1010,1030))
fig,ax=plt.subplots(1,1,figsize=(20,.3))
ax.plot(data.dtUTC,np.ones(np.shape(data.dtUTC)),'k.')
ax.set_xlim((dt.datetime(year,1,1),dt.datetime(year,12,31)))
ax.set_title('Data Timing')
ax.yaxis.set_visible(False)
def tsplot(ax,svar,tvar):
limsS=(0,36)
limsT=(5,20)
ss,tt=np.meshgrid(np.linspace(limsS[0],limsS[1],20),np.linspace(limsT[0],limsT[1],20))
rho=gsw.rho(ss,tt,np.zeros(np.shape(ss)))
r=ax.contour(ss,tt,rho,colors='k')
ps1=ax.plot(dJDF[svar],dJDF[tvar],'b.',label='SJDF')
ps2=ax.plot(dSJGI[svar],dSJGI[tvar],'c.',label='SJGI')
ps3=ax.plot(dSOG[svar],dSOG[tvar],'y.',label='SOG')
ps4=ax.plot(dNSOG[svar],dNSOG[tvar],'m.',label='NSOG')
l=ax.legend(handles=[ps1[0],ps2[0],ps3[0],ps4[0]],bbox_to_anchor=(1.55,1))
ax.set_ylim(limsT)
ax.set_xlim(limsS)
ax.set_ylabel('$\Theta$ ($^{\circ}$C)')
ax.set_xlabel('S$_A$ (g kg$^{-1}$)')
ax.set_aspect((limsS[1]-limsS[0])/(limsT[1]-limsT[0]))
return
fig,ax=plt.subplots(1,2,figsize=(16,4))
tsplot(ax[0],'AbsSal','ConsT')
ax[0].set_title('Observed')
tsplot(ax[1],'mod_vosaline','mod_votemper')
ax[1].set_title('Modelled')
Text(0.5, 1.0, 'Modelled')
# save stats dict to json file:
with open('vET-HC1905-DFO-NutChlPhys-stats.json', 'w') as fstat:
json.dump(statsDict, fstat, indent=4);
tbl,tdf=et.displayStats(statsDict[year],level='Variable',suborder=list(statsubs.keys()))
tbl
Bias | N | RMSE | WSS | |||
---|---|---|---|---|---|---|
Variable | Subset | |||||
Chl | 0 | z < 15 m | -1.54052 | 309 | 5.62303 | 0.41881 |
1 | 15 m < z < 22 m | -0.518307 | 116 | 1.78706 | 0.424334 | |
2 | z >= 22 m | -1.88707 | 3 | 2.57095 | 0.534625 | |
3 | z > 50 m | nan | 0 | nan | nan | |
4 | all | -1.2659 | 428 | 4.87229 | 0.452938 | |
5 | z < 15 m, JFM | -2.79809 | 40 | 7.79091 | 0.29813 | |
6 | z < 15 m, Apr | -5.43267 | 64 | 8.35052 | 0.553279 | |
7 | z < 15 m, MJJA | 0.521975 | 126 | 4.23393 | 0.36672 | |
8 | z < 15 m, SOND | -1.04019 | 79 | 2.80295 | 0.403642 | |
Chl log10 | 0 | z < 15 m | -0.116013 | 309 | 0.566884 | 0.514615 |
1 | 15 m < z < 22 m | -0.0796708 | 116 | 0.479648 | 0.543346 | |
2 | z >= 22 m | -0.441195 | 3 | 0.492631 | 0.611506 | |
3 | z > 50 m | nan | 0 | nan | nan | |
4 | all | -0.108443 | 428 | 0.544116 | 0.617788 | |
5 | z < 15 m, JFM | -0.309659 | 40 | 0.625575 | 0.386615 | |
6 | z < 15 m, Apr | -0.514023 | 64 | 0.695956 | 0.6122 | |
7 | z < 15 m, MJJA | 0.200752 | 126 | 0.475027 | 0.411375 | |
8 | z < 15 m, SOND | -0.200747 | 79 | 0.553642 | 0.408978 | |
Density | 0 | z < 15 m | -0.10353 | 655 | 1.27226 | 0.938817 |
1 | 15 m < z < 22 m | 0.0443576 | 158 | 0.454994 | 0.940907 | |
2 | z >= 22 m | 0.118125 | 1022 | 0.234376 | 0.993265 | |
3 | z > 50 m | 0.12618 | 647 | 0.209889 | 0.991867 | |
4 | all | 0.032654 | 1835 | 0.791324 | 0.97542 | |
5 | z < 15 m, JFM | -0.333153 | 54 | 1.18491 | 0.835796 | |
6 | z < 15 m, Apr | -0.100339 | 219 | 0.721694 | 0.933834 | |
7 | z < 15 m, MJJA | -0.0390939 | 274 | 1.68556 | 0.926471 | |
8 | z < 15 m, SOND | -0.158668 | 108 | 0.922291 | 0.917526 | |
NO3 | 0 | z < 15 m | 1.21132 | 508 | 6.09019 | 0.837433 |
1 | 15 m < z < 22 m | -0.211189 | 137 | 3.84811 | 0.840336 | |
2 | z >= 22 m | -0.739887 | 945 | 2.60958 | 0.893272 | |
3 | z > 50 m | -1.26612 | 579 | 2.20042 | 0.892088 | |
4 | all | -0.0709275 | 1590 | 4.1441 | 0.923395 | |
5 | z < 15 m, JFM | -2.65867 | 42 | 3.66444 | 0.822522 | |
6 | z < 15 m, Apr | 5.87382 | 185 | 7.35724 | 0.642506 | |
7 | z < 15 m, MJJA | -1.51093 | 184 | 4.70592 | 0.878824 | |
8 | z < 15 m, SOND | -0.841547 | 97 | 6.57179 | 0.786077 | |
Salinity | 0 | z < 15 m | -0.0904083 | 655 | 1.5314 | 0.931232 |
1 | 15 m < z < 22 m | 0.0580848 | 158 | 0.530179 | 0.94627 | |
2 | z >= 22 m | 0.132994 | 1022 | 0.28146 | 0.990731 | |
3 | z > 50 m | 0.147509 | 647 | 0.258548 | 0.989591 | |
4 | all | 0.0468009 | 1835 | 0.951543 | 0.967239 | |
5 | z < 15 m, JFM | -0.50658 | 54 | 1.58329 | 0.827336 | |
6 | z < 15 m, Apr | -0.13077 | 219 | 0.881454 | 0.935913 | |
7 | z < 15 m, MJJA | 0.0577092 | 274 | 2.00616 | 0.924753 | |
8 | z < 15 m, SOND | -0.176257 | 108 | 1.08788 | 0.891446 | |
Temperature | 0 | z < 15 m | 0.153269 | 655 | 1.10507 | 0.968584 |
1 | 15 m < z < 22 m | -0.0160375 | 158 | 0.511579 | 0.961782 | |
2 | z >= 22 m | -0.104113 | 1022 | 0.368898 | 0.96031 | |
3 | z > 50 m | -0.078967 | 647 | 0.354609 | 0.954128 | |
4 | all | -0.00465754 | 1835 | 0.730905 | 0.975418 | |
5 | z < 15 m, JFM | -0.520092 | 54 | 0.700092 | 0.549084 | |
6 | z < 15 m, Apr | -0.00468493 | 219 | 0.416628 | 0.842328 | |
7 | z < 15 m, MJJA | 0.429239 | 274 | 1.50409 | 0.911308 | |
8 | z < 15 m, SOND | 0.110097 | 108 | 1.03424 | 0.932098 | |
dSi | 0 | z < 15 m | 1.42277 | 508 | 13.4813 | 0.6459 |
1 | 15 m < z < 22 m | -3.897 | 137 | 10.3804 | 0.684622 | |
2 | z >= 22 m | -5.30834 | 944 | 8.63213 | 0.776742 | |
3 | z > 50 m | -5.8247 | 578 | 8.18832 | 0.823395 | |
4 | all | -3.03473 | 1589 | 10.567 | 0.781749 | |
5 | z < 15 m, JFM | -8.49306 | 42 | 9.46056 | 0.771748 | |
6 | z < 15 m, Apr | 14.7136 | 185 | 18.1407 | 0.478543 | |
7 | z < 15 m, MJJA | -6.8016 | 184 | 10.3377 | 0.691526 | |
8 | z < 15 m, SOND | -4.0315 | 97 | 9.09493 | 0.502836 |