In [1]:
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
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>''')
Out[2]:
In [3]:
PATH= '/results2/SalishSea/nowcast-green.201905/'
year=2007
In [4]:
# Parameters
year = 2015
In [5]:
display(Markdown('''# Year: '''+ str(year)))

Year: 2015

Yearly model-data comparisons of nutrients, chlorophyll, temperature and salinity between 201905 runs and DFO observations

Define date range and load observations

In [6]:
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()
569 data points
Out[6]:
Year Month Day Lat Lon Pressure Depth N Si Chlorophyll_Extracted ConsT AbsSal
0 2015.0 2.0 11.0 48.300833 -124.000333 1.9 None 15.31 32.14 NaN 9.859421 29.227507
1 2015.0 2.0 11.0 48.300833 -124.000333 6.6 None 17.13 33.90 2.57 9.777243 29.484341
2 2015.0 2.0 11.0 48.300833 -124.000333 6.7 None NaN NaN NaN 9.771987 29.484839
3 2015.0 2.0 11.0 48.300833 -124.000333 11.0 None NaN NaN NaN 9.439995 30.144549
4 2015.0 2.0 11.0 48.300833 -124.000333 11.0 None 20.62 37.65 NaN 9.433733 30.157913
In [7]:
data=et.matchData(df1,filemap,fdict,start_date,end_date,'nowcast',PATH,1,quiet=True);
In [8]:
# 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']))
In [9]:
# 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 01jan15 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.88
   zz_alpha_b_d =  0.0
In [10]:
# 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']
In [11]:
# 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()};
In [12]:
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'][:,:])
In [13]:
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
In [14]:
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)),:]
In [15]:
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})

Nitrate

In [16]:
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
Out[16]:
Bias N RMSE WSS
Subset
0 z < 15 m -0.69409 127 4.84458 0.889856
1 15 m < z < 22 m -0.49358 38 3.98725 0.6468
2 z >= 22 m -2.13697 368 2.80707 0.767835
3 z > 50 m -2.37822 286 2.82053 0.728406
4 all -1.67601 533 3.48799 0.934467
5 z < 15 m, JFM -1.477 8 3.0851 0.792953
6 z < 15 m, Apr -1.77523 67 5.52535 0.743014
7 z < 15 m, MJJA -1.12733 29 3.46787 0.896596
8 z < 15 m, SOND 3.2739 23 4.71031 0.840876
In [17]:
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');
In [18]:
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)
In [19]:
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))

Dissolved Silica

In [20]:
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
Out[20]:
Bias N RMSE WSS
Subset
0 z < 15 m -3.37252 127 8.9829 0.785152
1 15 m < z < 22 m -2.93607 38 8.0882 0.646524
2 z >= 22 m -6.38061 368 7.96496 0.788195
3 z > 50 m -6.87779 286 8.2254 0.790252
4 all -5.41828 533 8.22754 0.85954
5 z < 15 m, JFM -4.86894 8 6.22762 0.78727
6 z < 15 m, Apr -0.781434 67 9.99822 0.711514
7 z < 15 m, MJJA -6.54298 29 8.17587 0.792911
8 z < 15 m, SOND -6.40246 23 7.52265 0.598907
In [21]:
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');
In [22]:
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)
In [23]:
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))

Profiles of NO3 and Dissolved Silica

In [24]:
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')
Out[24]:
Text(0.5, 1.0, 'dSi')

dSi:NO3 Ratios

In [25]:
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-')
Out[25]:
<matplotlib.legend.Legend at 0x7fee75639220>
In [26]:
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()
Out[26]:
<matplotlib.legend.Legend at 0x7fee505e1910>

Chlorophyll

In [27]:
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
Out[27]:
Variable Chl Chl log10
Bias N RMSE WSS Bias N RMSE WSS
Subset
0 z < 15 m -0.114974 90 4.03009 0.521973 0.0107832 90 0.472198 0.537946
1 15 m < z < 22 m -0.38143 34 1.02139 0.632218 -0.0871032 34 0.295501 0.728333
2 z >= 22 m 0.0577223 1 0.0577223 0 0.0547041 1 0.0547041 0
3 z > 50 m nan 0 nan nan nan 0 nan nan
4 all -0.186068 125 3.46089 0.580621 -0.0154906 125 0.429319 0.676774
5 z < 15 m, JFM -0.859972 3 1.07978 0.552293 -0.34048 3 0.361261 0.561712
6 z < 15 m, Apr -0.187168 46 4.96997 0.508131 0.0641398 46 0.476528 0.553956
7 z < 15 m, MJJA 1.39006 23 2.9852 0.506031 0.284276 23 0.475082 0.490188
8 z < 15 m, SOND -1.72941 18 2.55016 0.436722 -0.416492 18 0.473701 0.465577
In [28]:
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');
In [29]:
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');

Conservative Temperature

In [30]:
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
Out[30]:
Bias N RMSE WSS
Subset
0 z < 15 m -0.132824 136 0.89454 0.977581
1 15 m < z < 22 m -0.159216 40 0.554603 0.934825
2 z >= 22 m 0.0553272 390 0.35997 0.933263
3 z > 50 m 0.0955125 303 0.3517 0.930767
4 all -0.00504425 566 0.550725 0.977906
5 z < 15 m, JFM -0.233373 13 0.493924 0.740153
6 z < 15 m, Apr -0.0808868 69 0.396238 0.748258
7 z < 15 m, MJJA 0.0112905 29 1.68229 0.904618
8 z < 15 m, SOND -0.391056 25 0.714121 0.922247
In [31]:
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');
In [32]:
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)

Reference Salinity

In [33]:
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
Out[33]:
Bias N RMSE WSS
Subset
0 z < 15 m -0.209793 136 1.28597 0.972119
1 15 m < z < 22 m 0.113384 40 0.478224 0.923512
2 z >= 22 m 0.0939602 390 0.262561 0.985888
3 z > 50 m 0.0901776 303 0.224463 0.988481
4 all 0.0223463 566 0.678989 0.983453
5 z < 15 m, JFM 0.320075 13 0.682321 0.383104
6 z < 15 m, Apr -0.226097 69 0.89296 0.95965
7 z < 15 m, MJJA -0.922775 29 2.21028 0.960429
8 z < 15 m, SOND 0.386733 25 0.941498 0.859428
In [34]:
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');
In [35]:
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)

Density

In [36]:
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
Out[36]:
Bias N RMSE WSS
Subset
0 z < 15 m -0.135476 136 1.08503 0.97379
1 15 m < z < 22 m 0.113226 40 0.437475 0.909251
2 z >= 22 m 0.0645063 390 0.224873 0.990031
3 z > 50 m 0.0552939 303 0.189724 0.991191
4 all 0.019897 566 0.575545 0.98673
5 z < 15 m, JFM 0.280876 13 0.577165 0.331427
6 z < 15 m, Apr -0.163554 69 0.712289 0.957209
7 z < 15 m, MJJA -0.691679 29 1.90624 0.958068
8 z < 15 m, SOND 0.37071 25 0.784719 0.882801
In [37]:
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');