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

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

In [3]:

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

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

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

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

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

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

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

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

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

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

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,
                      'SJDF':dJDF,
                      'SJGI':dSJGI,
                      'SOG':dSOG,
                      'NSOG':dNSOG})

Nitrate¶

In [14]:

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


		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
9	SJDF	-1.1742	64	3.93919	0.908224
10	SJGI	-2.17101	178	3.31845	0.853285
11	SOG	-1.51148	169	3.72604	0.949294
12	NSOG	-1.44494	122	3.11654	0.958259

In [15]:

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

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

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

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


		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
9	SJDF	0.10345	64	6.62626	0.813823
10	SJGI	-5.22237	178	6.66665	0.696064
11	SOG	-8.12679	169	9.69456	0.843495
12	NSOG	-4.84883	122	8.81295	0.899565

In [19]:

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

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

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

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

Text(0.5, 1.0, 'dSi')

dSi:NO3 Ratios¶

In [23]:

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

<matplotlib.legend.Legend at 0x7f08fcde4e50>

In [24]:

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

<matplotlib.legend.Legend at 0x7f08e2e02610>

Chlorophyll¶

In [25]:

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


	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
9	SJDF	-2.70618	12	3.85253	0.552701	-0.364343	12	0.481407	0.640331
10	SJGI	-0.479308	40	2.66698	0.606717	-0.0608658	40	0.392601	0.684754
11	SOG	0.279349	39	2.65934	0.574439	-0.0265022	39	0.407916	0.685107
12	NSOG	0.51451	34	4.72414	0.60584	0.173647	34	0.472896	0.694844

In [26]:

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

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

In [28]:

fig, ax = plt.subplots(1,1,figsize = (14,6))
obsvar='l10_obsChl'; modvar='l10_modChl'
with nc.Dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc') as fm:
    bathy=np.sum(fm.variables['tmask'][0,:,:,:]*fm.variables['e3t_0'][0,:,:,:],0)
data['bathy']=[bathy[j,i] for j,i in zip(data['j'],data['i'])]
ax.plot(data['bathy'],data['mod_Chl']-data['Chl'],'k.')
ax.set_title('Chl Residuals ($\mu$g/L) vs Model Depth');
ax.set_xlabel('Model Water Column Depth (m)')
ax.set_ylabel('Mod - Obs Chl ($\mu$g/L)');

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.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
9	SJDF	0.337086	77	0.504624	0.916498
10	SJGI	0.144091	188	0.473375	0.968285
11	SOG	-0.129842	176	0.678693	0.977637
12	NSOG	-0.264381	125	0.48051	0.985002

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.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
9	SJDF	0.0265695	77	0.470068	0.972096
10	SJGI	-0.0966657	188	0.514047	0.97831
11	SOG	0.0934705	176	0.793483	0.985704
12	NSOG	0.098596	125	0.817001	0.971554

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.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
9	SJDF	-0.0312604	77	0.402895	0.981746
10	SJGI	-0.0964276	188	0.452695	0.98222
11	SOG	0.0968652	176	0.687026	0.98829
12	NSOG	0.117991	125	0.653542	0.981222

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

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

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

Out[41]:


			Bias	N	RMSE	WSS
Variable		Subset
Chl	0	z < 15 m	-0.114974	90	4.03009	0.521973
	1	15 m < z < 22 m	-0.38143	34	1.02139	0.632218
	2	z >= 22 m	0.0577223	1	0.0577223	0
	3	z > 50 m	nan	0	nan	nan
	4	all	-0.186068	125	3.46089	0.580621
	5	z < 15 m, JFM	-0.859972	3	1.07978	0.552293
	6	z < 15 m, Apr	-0.187168	46	4.96997	0.508131
	7	z < 15 m, MJJA	1.39006	23	2.9852	0.506031
	8	z < 15 m, SOND	-1.72941	18	2.55016	0.436722
	9	SJDF	-2.70618	12	3.85253	0.552701
	10	SJGI	-0.479308	40	2.66698	0.606717
	11	SOG	0.279349	39	2.65934	0.574439
	12	NSOG	0.51451	34	4.72414	0.60584
Chl log10	0	z < 15 m	0.0107832	90	0.472198	0.537946
	1	15 m < z < 22 m	-0.0871032	34	0.295501	0.728333
	2	z >= 22 m	0.0547041	1	0.0547041	0
	3	z > 50 m	nan	0	nan	nan
	4	all	-0.0154906	125	0.429319	0.676774
	5	z < 15 m, JFM	-0.34048	3	0.361261	0.561712
	6	z < 15 m, Apr	0.0641398	46	0.476528	0.553956
	7	z < 15 m, MJJA	0.284276	23	0.475082	0.490188
	8	z < 15 m, SOND	-0.416492	18	0.473701	0.465577
	9	SJDF	-0.364343	12	0.481407	0.640331
	10	SJGI	-0.0608658	40	0.392601	0.684754
	11	SOG	-0.0265022	39	0.407916	0.685107
	12	NSOG	0.173647	34	0.472896	0.694844
Density	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
	9	SJDF	-0.0312604	77	0.402895	0.981746
	10	SJGI	-0.0964276	188	0.452695	0.98222
	11	SOG	0.0968652	176	0.687026	0.98829
	12	NSOG	0.117991	125	0.653542	0.981222
NO3	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
	9	SJDF	-1.1742	64	3.93919	0.908224
	10	SJGI	-2.17101	178	3.31845	0.853285
	11	SOG	-1.51148	169	3.72604	0.949294
	12	NSOG	-1.44494	122	3.11654	0.958259
Salinity	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
	9	SJDF	0.0265695	77	0.470068	0.972096
	10	SJGI	-0.0966657	188	0.514047	0.97831
	11	SOG	0.0934705	176	0.793483	0.985704
	12	NSOG	0.098596	125	0.817001	0.971554
Temperature	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
	9	SJDF	0.337086	77	0.504624	0.916498
	10	SJGI	0.144091	188	0.473375	0.968285
	11	SOG	-0.129842	176	0.678693	0.977637
	12	NSOG	-0.264381	125	0.48051	0.985002
dSi	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
	9	SJDF	0.10345	64	6.62626	0.813823
	10	SJGI	-5.22237	178	6.66665	0.696064
	11	SOG	-8.12679	169	9.69456	0.843495
	12	NSOG	-4.84883	122	8.81295	0.899565