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

PATH= '/results2/SalishSea/nowcast-green.201905/'
year=2015

In [4]:

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

start_date = dt.datetime(year,1,1)
end_date = dt.datetime(year,3,1)
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()

48 data points

Out[5]:

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

data=et.matchData(df1,filemap,fdict,start_date,end_date,'nowcast',PATH,1,quiet=True);

In [6]:

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

# report coefficients used in 
try:
    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'])
except:
    print('namelist_smelt_cfg not present in results file')

Parameter values from 01jan07 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 [9]:

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

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

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

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

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

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

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


		Bias	N	RMSE	WSS
	Subset
0	z < 15 m	-0.845001	230	4.8863	0.887092
1	15 m < z < 22 m	-0.287201	72	3.71001	0.558504
2	z >= 22 m	-0.682046	758	2.95439	0.744281
3	z > 50 m	-0.689607	589	2.7814	0.711584
4	all	-0.690585	1060	3.51528	0.894704
5	z < 15 m, JFM	-0.119254	3	2.54268	0
6	z < 15 m, Apr	2.68279	60	5.13338	0.530026
7	z < 15 m, MJJA	-2.1997	59	5.3411	0.873868
8	z < 15 m, SOND	-2.08498	108	4.5214	0.789183

In [16]:

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

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

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

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


		Bias	N	RMSE	WSS
	Subset
0	z < 15 m	-7.60506	230	14.4057	0.693065
1	15 m < z < 22 m	-5.77137	72	10.4514	0.483425
2	z >= 22 m	-5.54415	758	8.32201	0.725874
3	z > 50 m	-5.48	589	8.0562	0.731126
4	all	-6.00676	1060	10.0982	0.765622
5	z < 15 m, JFM	-1.26711	3	2.17827	0.374251
6	z < 15 m, Apr	5.89464	60	11.7713	0.277451
7	z < 15 m, MJJA	-14.5246	59	17.9863	0.621507
8	z < 15 m, SOND	-11.5008	108	13.7153	0.514492

In [20]:

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

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

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

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

Text(0.5, 1.0, 'dSi')

dSi:NO3 Ratios¶

In [24]:

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

<matplotlib.legend.Legend at 0x7f7284f95af0>

In [25]:

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

<matplotlib.legend.Legend at 0x7f72843c0520>

Chlorophyll¶

In [26]:

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


	Variable	Chl				Chl log10
		Bias	N	RMSE	WSS	Bias	N	RMSE	WSS
	Subset
0	z < 15 m	-1.09113	164	3.54545	0.489677	-0.189204	164	0.446323	0.73195
1	15 m < z < 22 m	-0.471902	58	1.07406	0.670933	-0.164835	58	0.363038	0.628456
2	z >= 22 m	-0.0490161	27	0.489433	0.324892	0.0152399	27	0.393419	0.469347
3	z > 50 m	0.0818484	4	0.0858041	0.379112	0.206992	4	0.212853	0.41876
4	all	-0.833894	249	2.92812	0.520093	-0.16136	249	0.422712	0.74179
5	z < 15 m, JFM	-5.56901	3	5.77952	0.295425	-0.874356	3	0.890711	0.160394
6	z < 15 m, Apr	-2.6436	47	5.98685	0.394372	-0.234814	47	0.395813	0.648459
7	z < 15 m, MJJA	0.0180185	40	2.17551	0.622982	0.00370198	40	0.275985	0.582558
8	z < 15 m, SOND	-0.523109	74	1.08681	0.467974	-0.236734	74	0.518303	0.431156

In [27]:

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

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

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


		Bias	N	RMSE	WSS
	Subset
0	z < 15 m	-0.0926801	236	0.706706	0.968923
1	15 m < z < 22 m	-0.151182	74	0.477925	0.927271
2	z >= 22 m	-0.146105	793	0.353202	0.94665
3	z > 50 m	-0.122045	619	0.315127	0.958586
4	all	-0.135015	1103	0.460298	0.964884
5	z < 15 m, JFM	-0.447043	3	0.675534	0.0961237
6	z < 15 m, Apr	-0.142572	61	0.196061	0.619533
7	z < 15 m, MJJA	0.010624	60	1.11232	0.936198
8	z < 15 m, SOND	-0.111356	112	0.596994	0.952276

In [30]:

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

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

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


		Bias	N	RMSE	WSS
	Subset
0	z < 15 m	-0.154293	236	1.54972	0.907624
1	15 m < z < 22 m	0.0786495	74	0.362279	0.973817
2	z >= 22 m	0.0998034	793	0.299041	0.986938
3	z > 50 m	0.0944484	619	0.296436	0.986564
4	all	0.0440172	1103	0.766132	0.960376
5	z < 15 m, JFM	-1.06164	3	1.74162	0.106754
6	z < 15 m, Apr	-0.0262992	61	0.735035	0.939549
7	z < 15 m, MJJA	-0.00339569	60	2.34858	0.903338
8	z < 15 m, SOND	-0.280538	112	1.31537	0.887514

In [33]:

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

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

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


		Bias	N	RMSE	WSS
	Subset
0	z < 15 m	-0.107181	236	1.23291	0.919497
1	15 m < z < 22 m	0.0831493	74	0.325284	0.967282
2	z >= 22 m	0.0993469	793	0.257235	0.989164
3	z > 50 m	0.0916767	619	0.249632	0.988553
4	all	0.054071	1103	0.616368	0.971349
5	z < 15 m, JFM	-0.774649	3	1.2803	0.10581
6	z < 15 m, Apr	-0.00113046	61	0.575546	0.938949
7	z < 15 m, MJJA	-8.30478e-05	60	1.89455	0.915124
8	z < 15 m, SOND	-0.204436	112	1.02754	0.892563

In [36]:

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');

In [37]:

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)

Temperature-Salinity by Region¶

In [38]:

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

In [39]:

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')

Out[39]:

Text(0.5, 1.0, 'Modelled')

In [41]:

# save stats dict to json file:
with open('vET-HC1905-DFO-NutChlPhys-stats.json', 'w') as fstat:
    json.dump(statsDict, fstat, indent=4);

Display All Stats¶

In [42]:

tbl,tdf=et.displayStats(statsDict[year],level='Variable',suborder=list(statsubs.keys()))
tbl

Out[42]:


			Bias	N	RMSE	WSS
Variable		Subset
Chl	0	z < 15 m	-1.09113	164	3.54545	0.489677
	1	15 m < z < 22 m	-0.471902	58	1.07406	0.670933
	2	z >= 22 m	-0.0490161	27	0.489433	0.324892
	3	z > 50 m	0.0818484	4	0.0858041	0.379112
	4	all	-0.833894	249	2.92812	0.520093
	5	z < 15 m, JFM	-5.56901	3	5.77952	0.295425
	6	z < 15 m, Apr	-2.6436	47	5.98685	0.394372
	7	z < 15 m, MJJA	0.0180185	40	2.17551	0.622982
	8	z < 15 m, SOND	-0.523109	74	1.08681	0.467974
Chl log10	0	z < 15 m	-0.189204	164	0.446323	0.73195
	1	15 m < z < 22 m	-0.164835	58	0.363038	0.628456
	2	z >= 22 m	0.0152399	27	0.393419	0.469347
	3	z > 50 m	0.206992	4	0.212853	0.41876
	4	all	-0.16136	249	0.422712	0.74179
	5	z < 15 m, JFM	-0.874356	3	0.890711	0.160394
	6	z < 15 m, Apr	-0.234814	47	0.395813	0.648459
	7	z < 15 m, MJJA	0.00370198	40	0.275985	0.582558
	8	z < 15 m, SOND	-0.236734	74	0.518303	0.431156
Density	0	z < 15 m	-0.107181	236	1.23291	0.919497
	1	15 m < z < 22 m	0.0831493	74	0.325284	0.967282
	2	z >= 22 m	0.0993469	793	0.257235	0.989164
	3	z > 50 m	0.0916767	619	0.249632	0.988553
	4	all	0.054071	1103	0.616368	0.971349
	5	z < 15 m, JFM	-0.774649	3	1.2803	0.10581
	6	z < 15 m, Apr	-0.00113046	61	0.575546	0.938949
	7	z < 15 m, MJJA	-8.30478e-05	60	1.89455	0.915124
	8	z < 15 m, SOND	-0.204436	112	1.02754	0.892563
NO3	0	z < 15 m	-0.845001	230	4.8863	0.887092
	1	15 m < z < 22 m	-0.287201	72	3.71001	0.558504
	2	z >= 22 m	-0.682046	758	2.95439	0.744281
	3	z > 50 m	-0.689607	589	2.7814	0.711584
	4	all	-0.690585	1060	3.51528	0.894704
	5	z < 15 m, JFM	-0.119254	3	2.54268	0
	6	z < 15 m, Apr	2.68279	60	5.13338	0.530026
	7	z < 15 m, MJJA	-2.1997	59	5.3411	0.873868
	8	z < 15 m, SOND	-2.08498	108	4.5214	0.789183
Salinity	0	z < 15 m	-0.154293	236	1.54972	0.907624
	1	15 m < z < 22 m	0.0786495	74	0.362279	0.973817
	2	z >= 22 m	0.0998034	793	0.299041	0.986938
	3	z > 50 m	0.0944484	619	0.296436	0.986564
	4	all	0.0440172	1103	0.766132	0.960376
	5	z < 15 m, JFM	-1.06164	3	1.74162	0.106754
	6	z < 15 m, Apr	-0.0262992	61	0.735035	0.939549
	7	z < 15 m, MJJA	-0.00339569	60	2.34858	0.903338
	8	z < 15 m, SOND	-0.280538	112	1.31537	0.887514
Temperature	0	z < 15 m	-0.0926801	236	0.706706	0.968923
	1	15 m < z < 22 m	-0.151182	74	0.477925	0.927271
	2	z >= 22 m	-0.146105	793	0.353202	0.94665
	3	z > 50 m	-0.122045	619	0.315127	0.958586
	4	all	-0.135015	1103	0.460298	0.964884
	5	z < 15 m, JFM	-0.447043	3	0.675534	0.0961237
	6	z < 15 m, Apr	-0.142572	61	0.196061	0.619533
	7	z < 15 m, MJJA	0.010624	60	1.11232	0.936198
	8	z < 15 m, SOND	-0.111356	112	0.596994	0.952276
dSi	0	z < 15 m	-7.60506	230	14.4057	0.693065
	1	15 m < z < 22 m	-5.77137	72	10.4514	0.483425
	2	z >= 22 m	-5.54415	758	8.32201	0.725874
	3	z > 50 m	-5.48	589	8.0562	0.731126
	4	all	-6.00676	1060	10.0982	0.765622
	5	z < 15 m, JFM	-1.26711	3	2.17827	0.374251
	6	z < 15 m, Apr	5.89464	60	11.7713	0.277451
	7	z < 15 m, MJJA	-14.5246	59	17.9863	0.621507
	8	z < 15 m, SOND	-11.5008	108	13.7153	0.514492