Oxygen, compare 202111 to 201905R¶

In [21]:

import cmocean.cm as cm
import datetime as dt
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import numpy as np
import pandas as pd
# import statfieldodels.api as sm
import xarray as xr

from salishsea_tools import viz_tools

Bot data from 2011 - 2019
PSF data from 2015 - 2017??

In [22]:

field = 'Dissolved Oxygen'
obsn = 'Oxygen_Dissolved'
modn = 'mod_dissolved_oxygen'
vmax = 420
vmin = 30
dbin = 10
scale = 260

cmap = 'pink_r'#cm.rain #cm.ice_r #'ocean_r' #'pink_r'

In [23]:

foramt = "{:.3}"
myformat = {'bias': foramt, 'rmse': foramt, 'crmse':foramt, 'swillmott': foramt, 'slopedev': foramt, 'const': foramt,
                         'systematic': foramt, 'nonsystematic':foramt, 
                    'spread': foramt}

def use_f_2(x):
    return "%.2f" % x

In [24]:

with xr.open_dataset('/home/sallen/MEOPAR/grid/mesh_mask201702.nc') as mesh:
    deptht = mesh.gdept_1d[0].values
bathy = xr.open_dataset('/home/sallen/MEOPAR/grid/bathymetry_202108.nc')

In [25]:

def bias(df, obs, mod):
    diffy = df[mod] - df[obs]
    return diffy.count(), diffy.mean()

In [26]:

def rmse(df, obs, mod):
    return (np.sqrt(((df[mod] - df[obs])**2).mean()))

In [27]:

def crmse(rmse, bias):
    return (np.sqrt(rmse**2 - bias**2))

In [28]:

def swillmott(df, obs, mod):
    meanobs = df[obs].mean()
    return (((df[mod] - df[obs])**2).sum()
            /(( (df[mod] - meanobs).abs() + (df[obs] - meanobs).abs() )**2).sum())

In [29]:

def slope_inter(df, obs, mod):
    X = df[obs]
    y = df[mod]
    X = sm.add_constant(X)

    # Fit and make the predictions by the model
    model = sm.OLS(y, X, missing='drop').fit()
    predictions = model.predict(X)
    
    nonsyst = np.sqrt(((y - predictions)**2).mean())
    systematic = np.sqrt(((predictions - df[obs])**2).mean())
    return model.params[obs], model.params['const'], systematic, nonsyst

In [30]:

def wolfram_perp(df, obsin, modin):
    
    mod = np.array(df[modin][(df[modin] == df[modin]) & (df[obsin] == df[obsin])])
    obs = np.array(df[obsin][(df[modin] == df[modin]) & (df[obsin] == df[obsin])])
    
    n = mod.shape[0]

    y2s = (mod**2).sum()
    ys = mod.sum()
    x2s = (obs**2).sum()
    xs = obs.sum()
    xys = (mod * obs).sum()
    B = 0.5 * ( (y2s - ys**2/n) - (x2s - xs**2/n)) / (xs * ys/n - xys)
    b1 = -B + np.sqrt(B*B + 1)
    a1 = (ys - b1*xs)/n
    predictions = a1 + b1 * df[obsin]
    
    nonsyst = np.sqrt(((df[modin] - predictions)**2).mean())
    systematic = np.sqrt(((predictions - df[obsin])**2).mean())

    return a1, b1, systematic, nonsyst

In [31]:

def spread(df, obs, mod):
    return 1 - ((df[mod] - df[mod].mean())**2).mean() / ((df[obs] - df[obs].mean())**2).mean() 

In [32]:

def read_pieces(pieces):
    temp1 = pd.read_csv(pieces[0])
    for piece in pieces[1:]:
        nextpiece = pd.read_csv(piece)
        if 'ferry' in piece:
            nextpiece['k'] = 0
        temp1 = pd.concat([temp1, nextpiece], ignore_index=True)
    return temp1

In [33]:

def plot_and_stats(temp1, name, idepth, jdepth, fig, ax, whichdepths, reverse, string, boxes=False, box=None, boxname=None):
    if boxes:
        corn = box
        select = temp1[(temp1.k >= idepth) & (temp1.k <= jdepth) &
                              (temp1.j >= corn[0]) & (temp1.j <= corn[1]) &
                              (temp1.i >= corn[2]) & (temp1.i <= corn[3]) &
                              (temp1[obsn] == temp1[obsn])]
        shift_text = 0.94
    else:
        select = temp1[(temp1.k >= idepth) & (temp1.k < jdepth) & (temp1[obsn] == temp1[obsn])]
        shift_text = 1
    
    if reverse:
        one = modn
        two = obsn
    else:
        one = obsn
        two = modn
    
    counts, xedges, yedges, color = ax.hist2d(select[one], 
                                              select[two], 
                                           bins=np.arange(vmin, vmax, dbin), norm=LogNorm(), cmap=cmap);
    fig.colorbar(color, ax=ax)
   
    number, tbias = bias(select, one, two)
    trmse = rmse(select, one, two)
    tcrmse = crmse(trmse, tbias)
    tswillmott = swillmott(select, one, two)
#    m, c, syst, nonsyst = slope_inter(select, one, two)
    a1, b1, syst, nonsyst = wolfram_perp(select, one, two)
    tspread = spread(select, one, two)
    ax.plot([vmin, vmax], [vmin, vmax], 'k-');
    xr = np.arange(vmin, vmax, 0.5)
#    ax.plot(xr, c + m * xr, 'r-');
    ax.plot(xr, a1 + b1 * xr, 'r-')
    sc = scale/12
    sh = 2*sc-1
    bot = scale
    top = bot + 2*sh
    ax.arrow(sc+vmin, bot+vmin, 0, sh-np.abs(tbias)/2, head_width=0.5*sc, head_length=0.2*sc, length_includes_head=True)
    ax.arrow(sc+vmin, top+vmin, 0, -sh+np.abs(tbias)/2, head_width=0.5*sc, head_length=0.2*sc, length_includes_head=True)
    ax.arrow(2*sc+vmin, bot+vmin, 0, sh-syst/2, head_width=0.5*sc, head_length=0.2*sc, length_includes_head=True)
    ax.arrow(2*sc+vmin, top+vmin, 0, -sh+syst/2, head_width=0.5*sc, head_length=0.2*sc, length_includes_head=True)
    ax.arrow(3*sc+vmin, bot+vmin, 0, sh-nonsyst/2, head_width=0.5*sc, head_length=0.2*sc, length_includes_head=True)
    ax.arrow(3*sc+vmin, top+vmin, 0, -sh+nonsyst/2, head_width=0.5*sc, head_length=0.2*sc, length_includes_head=True);
    Cp2 = {'number': number,
            'bias': tbias,
          'rmse': trmse,
           'crmse': tcrmse,
          'swillmott': tswillmott,
          'slopedev': 1-b1,
          'const': a1,
          'systematic': syst,
          'nonsystematic': nonsyst,
          'spread': tspread}
    ax.text(0.8*sc+vmin, 0.9*bot*shift_text+vmin, 'bias', rotation=90)
    ax.text(1.8*sc+vmin, 0.72*bot*shift_text+vmin-shift_text, 'systematic', rotation=90)
    ax.text(2.8*sc+vmin, 0.6*bot*shift_text+vmin-shift_text, 'non-systematic', rotation=90)
    
    ax.set_title(f'{name}, {field} {whichdepths}{string}');
    dCp2 = pd.DataFrame(data=Cp2, index=[name])

    ax.set_ylabel('Model (µmol/L)')
    ax.set_xlabel('Obs (µmol/L)')
    return dCp2, counts

In [34]:

def highlight_max_min(s):
    '''
    highlight the maximum in a Series yellow.
    '''

    is_max = abs(s) == abs(s).max()
    is_min = abs(s) == abs(s).min()
    color = []
    for v, v2 in zip(is_max, is_min):
        if v:
            color.append('red')
        elif v2:
            color.append('darkgreen')
        else:
            color.append('black')
    return ['color: %s' % color[i] for i in range(len(is_max))]

In [35]:

def plot_allyears(years, years_psf, years_pug, years_ctdfo, years_onc, years_ferry, idepth, jdepth, whichdepths, reverse=False):
    if reverse:
        string = '_reverse'
    else:
        string = ''
    fig, axs = plt.subplots(1, 2, figsize=(12, 5))
    startyear = years[0] #min(years[0], years_psf[0], years_pug[0])#, years_onc[0])
    endyear = years [-1] #max(years[-1], years_psf[-1], years_pug[-1])#, years_onc[-1])
    fig.suptitle(f'Year {startyear}-{endyear}')
    pieces1 = []
    for year in years:
        pieces1.append(
              f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_201905R_bot_{year}0101_{year}1231.csv')
    for year in years_psf:
        pieces1.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_201905R_psf_{year}0101_{year}1231.csv')
    for year in years_pug:
        pieces1.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_201905R_pugts_{year}0101_{year}1231.csv')
    for year in years_ctdfo:
        pieces1.append(
            f'/ocean/atall/MOAD/ObsModel/202111/ObsModel_201905R_ctd_from_dfo_{year}0101_{year}1231.csv')
    for year in years_onc:
        pieces1.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_201905R_onc_{year}0101_{year}1231.csv')
    for year in years_ferry:
        pieces1.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_201905R_ferry_{year}0101_{year}1231.csv')
    temp1 = read_pieces(pieces1)
    d201905, counts1 = plot_and_stats(temp1, '201905R', idepth, jdepth, fig, axs[0], whichdepths, reverse, string)
    
    pieces2 = []
    for year in years:
        pieces2.append(
              f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_202111_bot_{year}0101_{year}1231.csv')
    for year in years_psf:
        pieces2.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_202111_psf_{year}0101_{year}1231.csv')
    for year in years_pug:
        pieces2.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_202111_pugts_{year}0101_{year}1231.csv')
    for year in years_ctdfo:
        pieces1.append(
            f'/ocean/atall/MOAD/ObsModel/202111/ObsModel_202111_ctd_from_dfo_{year}0101_{year}1231.csv')
    for year in years_onc:
        pieces2.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_202111_onc_{year}0101_{year}1231.csv')
    for year in years_ferry:
        pieces2.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_202111_ferry_{year}0101_{year}1231.csv')
    temp2 = read_pieces(pieces2)
    d202111, counts2 = plot_and_stats(temp2, '202111', idepth, jdepth, fig, axs[1], whichdepths, reverse, string)

    alltogether = pd.concat([d201905, d202111], axis=0)
    
    fig.savefig(f'{whichdepths}_{field}{string}_201905_202111_for_allyears.png')
    
    f = open(f'./{whichdepths}_{field}{string}_201905_202111_for_allyears_table.tex', 'w')
    f.write(alltogether.to_latex(column_format='lcccccccccc', 
            formatters=["{:d}".format, use_f_2, use_f_2, use_f_2, use_f_2, use_f_2, use_f_2, use_f_2, use_f_2, use_f_2]))
    f.close()
    return alltogether, counts1, counts2

All Depths¶

In [36]:

idepth = 0
jdepth = 39

results, counts1, counts2 = plot_allyears([2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019]
                                          , [2015, 2016, 2017], []
                                          , [2011, 2012, 2013, 2014, 2015, 2016]
                                       , [] #   , [2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022]
                                       , [] #   , [2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022]
                                          , idepth, jdepth, 'alldepths')
results.style.format(myformat).apply(highlight_max_min)

Out[36]:

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	547613	91.2	1.16e+02	71.0	0.603	0.136	1.15e+02	91.4	66.3	0.0333
202111	10405	-3.24	31.6	31.5	0.0644	0.099	16.0	7.35	29.4	0.168

Top 5 m¶

In [37]:

idepth = 0
jdepth = 5

results, counts1, counts2 = plot_allyears([2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019]
                                            , [2015, 2016, 2017]
                                          , [], [2011, 2012, 2013, 2014, 2015, 2016]
                                       , [] 
                                       , [], 
                                          idepth, jdepth, 'top5m')
results.style.format(myformat).apply(highlight_max_min)

Out[37]:

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	538497	93.1	1.16e+02	69.9	0.607	0.385	1.6e+02	95.3	57.6	0.104
202111	1289	-7.9	44.6	43.9	0.186	0.321	80.5	21.4	34.3	0.426

Below 5 m¶

In [38]:

idepth = 5
jdepth = 39

results, counts1, counts2 = plot_allyears([2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019]
                                            , [2015, 2016, 2017]
                                          , [], [2011, 2012, 2013, 2014, 2015, 2016]
                                       , [] 
                                       , [], 
                                          idepth, jdepth, 'below5m')
results.style.format(myformat).apply(highlight_max_min)

Out[38]:

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	9116	-21.8	41.0	34.7	0.113	-0.116	-43.0	22.8	36.2	-0.204
202111	9116	-2.58	29.4	29.2	0.0703	0.086	13.1	5.68	27.6	0.145

Regionals¶

In [39]:

def stats_only(temp1, name, idepth, jdepth, corn, boxname):
    select = temp1[(temp1.k >= idepth) & (temp1.k <= jdepth) &
                              (temp1.j >= corn[0]) & (temp1.j <= corn[1]) &
                              (temp1.i >= corn[2]) & (temp1.i <= corn[3])]
    number, tbias = bias(select, obsn, modn)
    if number > 0:
        trmse = rmse(select, obsn, modn)
        tcrmse = crmse(trmse, tbias)
        tswillmott = swillmott(select, obsn, modn)
        a1, b1, syst, nonsyst = wolfram_perp(select, obsn, modn)
        tspread = spread(select, obsn, modn)
        Cp2 = {'number': number,
            'bias': tbias,
          'rmse': trmse,
          'crmse': tcrmse,
          'swillmott': tswillmott,
          'slopedev': 1-b1,
          'const': a1,
          'systematic': syst,
          'nonsystematic': nonsyst,
          'spread': tspread,
          'region': boxname}
        dCp2 = pd.DataFrame(data=Cp2, index=[name])
    else:
        dCp2 = 0
    return number, dCp2

In [40]:

def stats_allregions(years, years_psf, years_pug, years_ctdfo, years_onc, years_ferry, idepth, jdepth, whichdepths):
    pieces1 = []
    for year in years:
        pieces1.append(
              f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_201905R_bot_{year}0101_{year}1231.csv'
             )
    for year in years_psf:
        pieces1.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_201905R_psf_{year}0101_{year}1231.csv')
    for year in years_pug:
        pieces1.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_201905R_pugts_{year}0101_{year}1231.csv')
    for year in years_ctdfo:
        pieces1.append(
            f'/ocean/atall/MOAD/ObsModel/202111/ObsModel_201905R_ctd_from_dfo_{year}0101_{year}1231.csv')
    for year in years_onc:
        pieces1.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_201905R_onc_{year}0101_{year}1231.csv')    
    for year in years_ferry:
        pieces1.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_201905R_ferry_{year}0101_{year}1231.csv')    
    temp1 = read_pieces(pieces1)
    
    number, alltogether = stats_only(temp1, '201905R', idepth, jdepth, boxes[0], boxnames[0])
    for boxname, corn in zip(boxnames[1:], boxes[1:]):
        number, d201905 = stats_only(temp1, '201905R', idepth, jdepth, corn, boxname)
        if number > 0:
            alltogether = pd.concat([alltogether, d201905], axis=0)

    pieces2 = []
    for year in years:
        pieces2.append(
              f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_202111_bot_{year}0101_{year}1231.csv'
             )
    for year in years_psf:
        pieces2.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_202111_psf_{year}0101_{year}1231.csv')
    for year in years_pug:
        pieces2.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_202111_pugts_{year}0101_{year}1231.csv')
    for year in years_ctdfo:
        pieces1.append(
            f'/ocean/atall/MOAD/ObsModel/202111/ObsModel_202111_ctd_from_dfo_{year}0101_{year}1231.csv')
    for year in years_onc:
        pieces2.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_202111_onc_{year}0101_{year}1231.csv')
    for year in years_ferry:
        pieces2.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_202111_ferry_{year}0101_{year}1231.csv')
    temp2 = read_pieces(pieces2)
   
    for boxname, corn in zip(boxnames, boxes):
        number, d202111 = stats_only(temp2, '202111', idepth, jdepth, corn, boxname)
        if number > 0:
            alltogether = pd.concat([alltogether, d202111], axis=0)
    f = open(f'./{whichdepths}_{field}_201905_202111_for_allregions_table.tex', 'w')
    f.write(alltogether.to_latex(column_format='lccccccccccl', 
            formatters=["{:d}".format, use_f_2, use_f_2, use_f_2, use_f_2, use_f_2, use_f_2, use_f_2, use_f_2, use_f_2, str.upper]))
    f.close()
    return number, alltogether

In [41]:

def stats_and_plot(years, years_psf, years_pug, years_ctdfo, years_onc, years_ferry, idepth, jdepth, box, boxname, whichdepths):
    fig, axs = plt.subplots(1, 4, figsize=(20, 4))
    fig.suptitle(f'{boxname} Year {years[0]}-{years[-1]}')
    
    pieces1 = []
    for year in years:
        pieces1.append(
              f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_201905R_bot_{year}0101_{year}1231.csv'
             )
    for year in years_psf:
        pieces1.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_201905R_psf_{year}0101_{year}1231.csv')
    for year in years_pug:
        pieces1.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_201905R_pugts_{year}0101_{year}1231.csv')
    for year in years_ctdfo:
        pieces1.append(
            f'/ocean/atall/MOAD/ObsModel/202111/ObsModel_201905R_ctd_from_dfo_{year}0101_{year}1231.csv')
    for year in years_onc:
        pieces1.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_201905R_onc_{year}0101_{year}1231.csv')   
    for year in years_ferry:
        pieces1.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_201905R_ferry_{year}0101_{year}1231.csv')   
    temp1 = read_pieces(pieces1)
    
    alltogether, counts = plot_and_stats(temp1, '201905R', idepth, jdepth, fig, axs[0], whichdepths, False, '', boxes=True, box=box, boxname=boxname)

    pieces2 = []
    for year in years:
        pieces2.append(
              f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_202111_bot_{year}0101_{year}1231.csv'
             )
    for year in years_psf:
        pieces2.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_202111_psf_{year}0101_{year}1231.csv')
    for year in years_pug:
        pieces2.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_202111_pugts_{year}0101_{year}1231.csv')
    for year in years_ctdfo:
        pieces1.append(
            f'/ocean/atall/MOAD/ObsModel/202111/ObsModel_202111_ctd_from_dfo_{year}0101_{year}1231.csv')
    for year in years_onc:
        pieces2.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_202111_onc_{year}0101_{year}1231.csv')
    for year in years_ferry:
        pieces2.append(
            f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_202111_ferry_{year}0101_{year}1231.csv')
    temp2 = read_pieces(pieces2)
    d202111, counts = plot_and_stats(temp2, '202111', idepth, jdepth, fig, axs[1], whichdepths, False, '', boxes=True, box=box, boxname=boxname)
    
    plot_profiles(temp1, temp2, axs[2], axs[3], box, boxname)
    
    
    alltogether = pd.concat([alltogether, d202111], axis=0)

In [42]:

def plot_profiles(temp1, temp2, ax, ax2, box, boxname):

    corn = box
    temp1['diff'] = temp1[modn] - temp1[obsn]
    select = temp1[(temp1.k >= idepth) & (temp1.k <= jdepth) &
                              (temp1.j >= corn[0]) & (temp1.j <= corn[1]) &
                              (temp1.i >= corn[2]) & (temp1.i <= corn[3]) &
                              (temp1[obsn] == temp1[obsn])]
    ks = (select.groupby('k')[obsn].mean().index[:])
    ax.plot(select.groupby('k')[obsn].mean(), deptht[ks], 'o-', c='tab:green', label='Observations')
    ax.plot(select.groupby('k')[modn].mean(), deptht[ks], 'o-', c='tab:blue', label='201905')
    ax2.plot(select.groupby('k')['diff'].mean(), deptht[ks], '-', c='tab:blue', label='201905')
    ax.fill_betweenx(deptht[ks], select[['k', modn]].groupby('k').quantile(q=0.25)[modn], select[['k', modn]].groupby('k').quantile(q=0.75)[modn], alpha=0.2, color='tab:blue')
    ax2.fill_betweenx(deptht[ks], select[['k', 'diff']].groupby('k').quantile(q=0.25)['diff'], select[['k', 'diff']].groupby('k').quantile(q=0.75)['diff'], alpha=0.2, color='tab:blue')
    
    temp2['diff'] = temp2[modn] - temp2[obsn]
    select = temp2[(temp2.k >= idepth) & (temp2.k <= jdepth) &
                              (temp2.j >= corn[0]) & (temp2.j <= corn[1]) &
                              (temp2.i >= corn[2]) & (temp2.i <= corn[3]) &
                              (temp2[obsn] == temp2[obsn])]
    
    ax.plot(select.groupby('k')[modn].mean(), deptht[ks], 'o-', c='tab:orange', label='202111')
    ax2.plot(select.groupby('k')['diff'].mean(), deptht[ks], '-', c='tab:orange', label='202111')
    ax.fill_betweenx(deptht[ks], select[['k', modn]].groupby('k').quantile(q=0.25)[modn], select[['k', modn]].groupby('k').quantile(q=0.75)[modn], alpha=0.2, color='tab:orange')
    ax.fill_betweenx(deptht[ks], select[['k', obsn]].groupby('k').quantile(q=0.25)[obsn], select[['k', obsn]].groupby('k').quantile(q=0.75)[obsn], alpha=0.2, color='tab:orange')
    ax2.fill_betweenx(deptht[ks], select[['k', 'diff']].groupby('k').quantile(q=0.25)['diff'], select[['k', 'diff']].groupby('k').quantile(q=0.75)['diff'], alpha=0.2, color='tab:orange')
    
    ax.invert_yaxis()
    ax2.invert_yaxis()
    
    ax.legend()
    ax2.legend()
    ax.set_title(f'{field} Profiles')
    ax2.set_title(f'{field} Difference Profiles')
    ax2.grid();
    ax.set_ylabel('Depth (m)')
    ax2.set_ylabel('Depth (m)')
    ax.set_xlabel('DO (µmol/L)')
    ax2.set_xlabel('DO (µmol/L)')

In [43]:

def plot_box(ax, corn, colour):
    ax.plot([corn[2], corn[3], corn[3], corn[2], corn[2]], 
            [corn[0], corn[0], corn[1], corn[1], corn[0]], '-', color=colour)

In [44]:

fig, ax = plt.subplots(1, 1, figsize=(5, 9))
mycmap = cm.deep
mycmap.set_bad('grey')
ax.pcolormesh(bathy['Bathymetry'], cmap=mycmap)
viz_tools.set_aspect(ax);
#Saanich_Inlet = [330, 390, 180, 220]
#plot_box(ax, Saanich_Inlet, 'm')
#Hood_Canal = [68, 150, 95, 200]
#plot_box(ax, Hood_Canal, 'm')
SoG_center = [450, 550, 200, 300]
plot_box(ax, SoG_center, 'b')
SoG_north = [650, 730, 100, 200]
plot_box(ax, SoG_north, 'g')
SoG_south = [320, 380, 280, 350]
plot_box(ax, SoG_south, 'k')
Haro_Boundary = [290, 350, 210, 280]
plot_box(ax, Haro_Boundary, 'g')
JdF_east = [200, 290, 150, 260]
plot_box(ax, JdF_east, 'b')
JdF_west = [250, 425, 25, 125]
plot_box(ax, JdF_west, 'k')
#PS_main = [20, 150, 200, 280]
PS_main = [0, 150, 100, 280]
plot_box(ax, PS_main, 'k')
SS_all = [0, 898, 0, 398]
boxes = [SS_all, SoG_south,SoG_center, SoG_north, PS_main, Haro_Boundary, JdF_west, JdF_east]
boxnames = ['SS_all','SoG_south', 'SoG_center', 'SoG_north', 'PS_main', 'Haro_Boundary',
           'JdF_west', 'JdF_east']
#boxnames = ['SS_all', 'Hood_Canal','Saanich_Inlet','SoG_south', 'SoG_center', 'SoG_north', 'PS_main', 'Haro_Boundary',
#           'JdF_west', 'JdF_east']
Fraser_plume = [380, 460, 260, 330]
plot_box(ax, Fraser_plume, 'm')

All Depths¶

In [46]:

idepth = 0
jdepth = 39

number, results = stats_allregions([2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019], 
                           [2015, 2016, 2017], [], [2011, 2012, 2013, 2014, 2015, 2016], [], 
                           [],
                           idepth, jdepth, 'alldepths')
for boxname in boxnames:
    if number > 0:
        print(f'\033[1m{boxname}\033[1m')
        display(results[results.region == boxname].loc[:, results.columns != 'region'].style.format(myformat).apply(highlight_max_min))

SS_all

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	547613	91.2	1.16e+02	71.0	0.603	0.136	1.15e+02	91.4	66.3	0.0989
202111	10405	-3.24	31.6	31.5	0.0644	0.099	16.0	7.35	29.4	0.0625

SoG_south

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	25184	62.8	80.8	50.8	0.486	-0.522	-42.8	66.7	63.2	-0.446
202111	700	1.76	17.9	17.8	0.0414	0.00419	2.61	1.77	17.8	-0.146

SoG_center

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	91380	1.18e+02	1.29e+02	52.5	0.635	0.274	1.65e+02	1.19e+02	45.5	0.379
202111	1680	-0.828	27.5	27.5	0.0502	0.1	17.9	6.6	25.5	0.00402

SoG_north

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	47261	97.9	1.15e+02	60.1	0.624	0.224	1.36e+02	98.5	53.6	0.0888
202111	954	0.0994	40.5	40.5	0.144	0.238	45.0	14.6	34.1	0.0801

PS_main

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	3119	64.0	89.5	62.5	0.638	-5.03	-9.06e+02	1.48e+02	1.5e+02	-4.67

Haro_Boundary

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	32612	34.5	55.4	43.3	0.357	-0.289	-18.9	36.4	49.1	-0.294
202111	691	1.99	28.5	28.4	0.0998	0.0751	16.2	4.13	27.2	0.119

JdF_west

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	23121	85.3	1.09e+02	67.9	0.533	0.601	1.82e+02	92.5	48.6	0.346
202111	912	-2.5	22.0	21.9	0.0249	-0.0506	-10.8	4.29	22.2	-0.0887

JdF_east

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	20388	35.6	57.8	45.5	0.325	0.158	64.8	36.5	41.7	0.168
202111	726	-2.86	18.0	17.8	0.0316	0.0565	7.68	4.13	17.1	0.128

In [47]:

stats_and_plot([2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019], [2015, 2016, 2017], 
                [], [2011, 2012, 2013, 2014, 2015, 2016], 
               [], 
               [],
               idepth, jdepth, SoG_south, 'SoG_south', 'alldepths')

In [48]:

stats_and_plot([2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019], [2015, 2016, 2017], 
                [], [2011, 2012, 2013, 2014, 2015, 2016], [], 
               [],
               idepth, jdepth, SoG_center, 'SoG_center', 'alldepths')

In [49]:

stats_and_plot([2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019], [2015, 2016, 2017], 
               [], [2011, 2012, 2013, 2014, 2015, 2016], [], [],
               idepth, jdepth, SoG_north, 'SoG_north', 'alldepths')

In [50]:

idepth = 0
jdepth = 39
stats_and_plot([2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019], [2015, 2016, 2017], 
               [], [2011, 2012, 2013, 2014, 2015, 2016], [], [],
               idepth, jdepth, PS_main, 'PS_main', 'alldepths')

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[50], line 3
      1 idepth = 0
      2 jdepth = 39
----> 3 stats_and_plot([2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019], [2015, 2016, 2017], 
      4                [], [2011, 2012, 2013, 2014, 2015, 2016], [], [],
      5                idepth, jdepth, PS_main, 'PS_main', 'alldepths')

Cell In[41], line 50, in stats_and_plot(years, years_psf, years_pug, years_ctdfo, years_onc, years_ferry, idepth, jdepth, box, boxname, whichdepths)
     47     pieces2.append(
     48         f'/data/sallen/results/MEOPAR/202111/ObsModel/ObsModel_202111_ferry_{year}0101_{year}1231.csv')
     49 temp2 = read_pieces(pieces2)
---> 50 d202111, counts = plot_and_stats(temp2, '202111', idepth, jdepth, fig, axs[1], whichdepths, False, '', boxes=True, box=box, boxname=boxname)
     52 plot_profiles(temp1, temp2, axs[2], axs[3], box, boxname)
     55 alltogether = pd.concat([alltogether, d202111], axis=0)

Cell In[33], line 23, in plot_and_stats(temp1, name, idepth, jdepth, fig, ax, whichdepths, reverse, string, boxes, box, boxname)
     18     two = modn
     20 counts, xedges, yedges, color = ax.hist2d(select[one], 
     21                                           select[two], 
     22                                        bins=np.arange(vmin, vmax, dbin), norm=LogNorm(), cmap=cmap);
---> 23 fig.colorbar(color, ax=ax)
     25 number, tbias = bias(select, one, two)
     26 trmse = rmse(select, one, two)

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/figure.py:1293, in FigureBase.colorbar(self, mappable, cax, ax, use_gridspec, **kwargs)
   1289 NON_COLORBAR_KEYS = ['fraction', 'pad', 'shrink', 'aspect', 'anchor',
   1290                      'panchor']
   1291 cb_kw = {k: v for k, v in kwargs.items() if k not in NON_COLORBAR_KEYS}
-> 1293 cb = cbar.Colorbar(cax, mappable, **cb_kw)
   1295 if not userax:
   1296     self.sca(current_ax)

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/_api/deprecation.py:384, in delete_parameter.<locals>.wrapper(*inner_args, **inner_kwargs)
    379 @functools.wraps(func)
    380 def wrapper(*inner_args, **inner_kwargs):
    381     if len(inner_args) <= name_idx and name not in inner_kwargs:
    382         # Early return in the simple, non-deprecated case (much faster than
    383         # calling bind()).
--> 384         return func(*inner_args, **inner_kwargs)
    385     arguments = signature.bind(*inner_args, **inner_kwargs).arguments
    386     if is_varargs and arguments.get(name):

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/colorbar.py:402, in Colorbar.__init__(self, ax, mappable, cmap, norm, alpha, values, boundaries, orientation, ticklocation, extend, spacing, ticks, format, drawedges, filled, extendfrac, extendrect, label, location)
    399 self.ticklocation = ticklocation
    401 self.set_label(label)
--> 402 self._reset_locator_formatter_scale()
    404 if np.iterable(ticks):
    405     self._locator = ticker.FixedLocator(ticks, nbins=len(ticks))

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/colorbar.py:1173, in Colorbar._reset_locator_formatter_scale(self)
   1167 def _reset_locator_formatter_scale(self):
   1168     """
   1169     Reset the locator et al to defaults.  Any user-hardcoded changes
   1170     need to be re-entered if this gets called (either at init, or when
   1171     the mappable normal gets changed: Colorbar.update_normal)
   1172     """
-> 1173     self._process_values()
   1174     self._locator = None
   1175     self._minorlocator = None

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/colorbar.py:1107, in Colorbar._process_values(self)
   1105     self.norm.vmin = 0
   1106     self.norm.vmax = 1
-> 1107 self.norm.vmin, self.norm.vmax = mtransforms.nonsingular(
   1108     self.norm.vmin, self.norm.vmax, expander=0.1)
   1109 if (not isinstance(self.norm, colors.BoundaryNorm) and
   1110         (self.boundaries is None)):
   1111     b = self.norm.inverse(b)

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/colors.py:1250, in Normalize.vmin(self, value)
   1248 if value != self._vmin:
   1249     self._vmin = value
-> 1250     self._changed()

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/colors.py:1278, in Normalize._changed(self)
   1273 def _changed(self):
   1274     """
   1275     Call this whenever the norm is changed to notify all the
   1276     callback listeners to the 'changed' signal.
   1277     """
-> 1278     self.callbacks.process('changed')

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/cbook/__init__.py:309, in CallbackRegistry.process(self, s, *args, **kwargs)
    307 except Exception as exc:
    308     if self.exception_handler is not None:
--> 309         self.exception_handler(exc)
    310     else:
    311         raise

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/cbook/__init__.py:96, in _exception_printer(exc)
     94 def _exception_printer(exc):
     95     if _get_running_interactive_framework() in ["headless", None]:
---> 96         raise exc
     97     else:
     98         traceback.print_exc()

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/cbook/__init__.py:304, in CallbackRegistry.process(self, s, *args, **kwargs)
    302 if func is not None:
    303     try:
--> 304         func(*args, **kwargs)
    305     # this does not capture KeyboardInterrupt, SystemExit,
    306     # and GeneratorExit
    307     except Exception as exc:

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/cm.py:665, in ScalarMappable.changed(self)
    660 def changed(self):
    661     """
    662     Call this whenever the mappable is changed to notify all the
    663     callbackSM listeners to the 'changed' signal.
    664     """
--> 665     self.callbacks.process('changed', self)
    666     self.stale = True

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/cbook/__init__.py:309, in CallbackRegistry.process(self, s, *args, **kwargs)
    307 except Exception as exc:
    308     if self.exception_handler is not None:
--> 309         self.exception_handler(exc)
    310     else:
    311         raise

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/cbook/__init__.py:96, in _exception_printer(exc)
     94 def _exception_printer(exc):
     95     if _get_running_interactive_framework() in ["headless", None]:
---> 96         raise exc
     97     else:
     98         traceback.print_exc()

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/cbook/__init__.py:304, in CallbackRegistry.process(self, s, *args, **kwargs)
    302 if func is not None:
    303     try:
--> 304         func(*args, **kwargs)
    305     # this does not capture KeyboardInterrupt, SystemExit,
    306     # and GeneratorExit
    307     except Exception as exc:

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/colorbar.py:514, in Colorbar.update_normal(self, mappable)
    511     self.norm = mappable.norm
    512     self._reset_locator_formatter_scale()
--> 514 self._draw_all()
    515 if isinstance(self.mappable, contour.ContourSet):
    516     CS = self.mappable

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/colorbar.py:549, in Colorbar._draw_all(self)
    543 self._short_axis().set_ticks([], minor=True)
    545 # Set self._boundaries and self._values, including extensions.
    546 # self._boundaries are the edges of each square of color, and
    547 # self._values are the value to map into the norm to get the
    548 # color:
--> 549 self._process_values()
    550 # Set self.vmin and self.vmax to first and last boundary, excluding
    551 # extensions:
    552 self.vmin, self.vmax = self._boundaries[self._inside][[0, -1]]

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/colorbar.py:1111, in Colorbar._process_values(self)
   1107 self.norm.vmin, self.norm.vmax = mtransforms.nonsingular(
   1108     self.norm.vmin, self.norm.vmax, expander=0.1)
   1109 if (not isinstance(self.norm, colors.BoundaryNorm) and
   1110         (self.boundaries is None)):
-> 1111     b = self.norm.inverse(b)
   1113 self._boundaries = np.asarray(b, dtype=float)
   1114 self._values = 0.5 * (self._boundaries[:-1] + self._boundaries[1:])

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/colors.py:1724, in _make_norm_from_scale.<locals>.Norm.inverse(self, value)
   1722 t_vmin, t_vmax = self._trf.transform([self.vmin, self.vmax])
   1723 if not np.isfinite([t_vmin, t_vmax]).all():
-> 1724     raise ValueError("Invalid vmin or vmax")
   1725 value, is_scalar = self.process_value(value)
   1726 rescaled = value * (t_vmax - t_vmin)

ValueError: Invalid vmin or vmax

Error in callback <function _draw_all_if_interactive at 0x7f4944cb9760> (for post_execute):

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/pyplot.py:120, in _draw_all_if_interactive()
    118 def _draw_all_if_interactive():
    119     if matplotlib.is_interactive():
--> 120         draw_all()

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/_pylab_helpers.py:132, in Gcf.draw_all(cls, force)
    130 for manager in cls.get_all_fig_managers():
    131     if force or manager.canvas.figure.stale:
--> 132         manager.canvas.draw_idle()

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/backend_bases.py:2082, in FigureCanvasBase.draw_idle(self, *args, **kwargs)
   2080 if not self._is_idle_drawing:
   2081     with self._idle_draw_cntx():
-> 2082         self.draw(*args, **kwargs)

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:400, in FigureCanvasAgg.draw(self)
    396 # Acquire a lock on the shared font cache.
    397 with RendererAgg.lock, \
    398      (self.toolbar._wait_cursor_for_draw_cm() if self.toolbar
    399       else nullcontext()):
--> 400     self.figure.draw(self.renderer)
    401     # A GUI class may be need to update a window using this draw, so
    402     # don't forget to call the superclass.
    403     super().draw()

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/artist.py:95, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     93 @wraps(draw)
     94 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 95     result = draw(artist, renderer, *args, **kwargs)
     96     if renderer._rasterizing:
     97         renderer.stop_rasterizing()

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     69     if artist.get_agg_filter() is not None:
     70         renderer.start_filter()
---> 72     return draw(artist, renderer)
     73 finally:
     74     if artist.get_agg_filter() is not None:

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/figure.py:3140, in Figure.draw(self, renderer)
   3137         # ValueError can occur when resizing a window.
   3139 self.patch.draw(renderer)
-> 3140 mimage._draw_list_compositing_images(
   3141     renderer, self, artists, self.suppressComposite)
   3143 for sfig in self.subfigs:
   3144     sfig.draw(renderer)

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/image.py:131, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    129 if not_composite or not has_images:
    130     for a in artists:
--> 131         a.draw(renderer)
    132 else:
    133     # Composite any adjacent images together
    134     image_group = []

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     69     if artist.get_agg_filter() is not None:
     70         renderer.start_filter()
---> 72     return draw(artist, renderer)
     73 finally:
     74     if artist.get_agg_filter() is not None:

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/axes/_base.py:3064, in _AxesBase.draw(self, renderer)
   3061 if artists_rasterized:
   3062     _draw_rasterized(self.figure, artists_rasterized, renderer)
-> 3064 mimage._draw_list_compositing_images(
   3065     renderer, self, artists, self.figure.suppressComposite)
   3067 renderer.close_group('axes')
   3068 self.stale = False

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/image.py:131, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    129 if not_composite or not has_images:
    130     for a in artists:
--> 131         a.draw(renderer)
    132 else:
    133     # Composite any adjacent images together
    134     image_group = []

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     69     if artist.get_agg_filter() is not None:
     70         renderer.start_filter()
---> 72     return draw(artist, renderer)
     73 finally:
     74     if artist.get_agg_filter() is not None:

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/collections.py:2074, in QuadMesh.draw(self, renderer)
   2071     ys = self.convert_yunits(offsets[:, 1])
   2072     offsets = np.column_stack([xs, ys])
-> 2074 self.update_scalarmappable()
   2076 if not transform.is_affine:
   2077     coordinates = self._coordinates.reshape((-1, 2))

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/collections.py:887, in Collection.update_scalarmappable(self)
    885         # pcolormesh, scatter, maybe others flatten their _A
    886         self._alpha = self._alpha.reshape(self._A.shape)
--> 887     self._mapped_colors = self.to_rgba(self._A, self._alpha)
    889 if self._face_is_mapped:
    890     self._facecolors = self._mapped_colors

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/cm.py:493, in ScalarMappable.to_rgba(self, x, alpha, bytes, norm)
    491 x = ma.asarray(x)
    492 if norm:
--> 493     x = self.norm(x)
    494 rgba = self.cmap(x, alpha=alpha, bytes=bytes)
    495 return rgba

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/colors.py:1711, in _make_norm_from_scale.<locals>.Norm.__call__(self, value, clip)
   1709 t_vmin, t_vmax = self._trf.transform([self.vmin, self.vmax])
   1710 if not np.isfinite([t_vmin, t_vmax]).all():
-> 1711     raise ValueError("Invalid vmin or vmax")
   1712 t_value -= t_vmin
   1713 t_value /= (t_vmax - t_vmin)

ValueError: Invalid vmin or vmax

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/IPython/core/formatters.py:340, in BaseFormatter.__call__(self, obj)
    338     pass
    339 else:
--> 340     return printer(obj)
    341 # Finally look for special method names
    342 method = get_real_method(obj, self.print_method)

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/IPython/core/pylabtools.py:152, in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
    149     from matplotlib.backend_bases import FigureCanvasBase
    150     FigureCanvasBase(fig)
--> 152 fig.canvas.print_figure(bytes_io, **kw)
    153 data = bytes_io.getvalue()
    154 if fmt == 'svg':

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/backend_bases.py:2342, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2336     renderer = _get_renderer(
   2337         self.figure,
   2338         functools.partial(
   2339             print_method, orientation=orientation)
   2340     )
   2341     with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2342         self.figure.draw(renderer)
   2344 if bbox_inches:
   2345     if bbox_inches == "tight":

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/artist.py:95, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     93 @wraps(draw)
     94 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 95     result = draw(artist, renderer, *args, **kwargs)
     96     if renderer._rasterizing:
     97         renderer.stop_rasterizing()

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     69     if artist.get_agg_filter() is not None:
     70         renderer.start_filter()
---> 72     return draw(artist, renderer)
     73 finally:
     74     if artist.get_agg_filter() is not None:

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/figure.py:3140, in Figure.draw(self, renderer)
   3137         # ValueError can occur when resizing a window.
   3139 self.patch.draw(renderer)
-> 3140 mimage._draw_list_compositing_images(
   3141     renderer, self, artists, self.suppressComposite)
   3143 for sfig in self.subfigs:
   3144     sfig.draw(renderer)

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/image.py:131, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    129 if not_composite or not has_images:
    130     for a in artists:
--> 131         a.draw(renderer)
    132 else:
    133     # Composite any adjacent images together
    134     image_group = []

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     69     if artist.get_agg_filter() is not None:
     70         renderer.start_filter()
---> 72     return draw(artist, renderer)
     73 finally:
     74     if artist.get_agg_filter() is not None:

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/axes/_base.py:3064, in _AxesBase.draw(self, renderer)
   3061 if artists_rasterized:
   3062     _draw_rasterized(self.figure, artists_rasterized, renderer)
-> 3064 mimage._draw_list_compositing_images(
   3065     renderer, self, artists, self.figure.suppressComposite)
   3067 renderer.close_group('axes')
   3068 self.stale = False

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/image.py:131, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    129 if not_composite or not has_images:
    130     for a in artists:
--> 131         a.draw(renderer)
    132 else:
    133     # Composite any adjacent images together
    134     image_group = []

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     69     if artist.get_agg_filter() is not None:
     70         renderer.start_filter()
---> 72     return draw(artist, renderer)
     73 finally:
     74     if artist.get_agg_filter() is not None:

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/collections.py:2074, in QuadMesh.draw(self, renderer)
   2071     ys = self.convert_yunits(offsets[:, 1])
   2072     offsets = np.column_stack([xs, ys])
-> 2074 self.update_scalarmappable()
   2076 if not transform.is_affine:
   2077     coordinates = self._coordinates.reshape((-1, 2))

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/collections.py:887, in Collection.update_scalarmappable(self)
    885         # pcolormesh, scatter, maybe others flatten their _A
    886         self._alpha = self._alpha.reshape(self._A.shape)
--> 887     self._mapped_colors = self.to_rgba(self._A, self._alpha)
    889 if self._face_is_mapped:
    890     self._facecolors = self._mapped_colors

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/cm.py:493, in ScalarMappable.to_rgba(self, x, alpha, bytes, norm)
    491 x = ma.asarray(x)
    492 if norm:
--> 493     x = self.norm(x)
    494 rgba = self.cmap(x, alpha=alpha, bytes=bytes)
    495 return rgba

File ~/conda_envs/analysis-abdoul/lib/python3.11/site-packages/matplotlib/colors.py:1711, in _make_norm_from_scale.<locals>.Norm.__call__(self, value, clip)
   1709 t_vmin, t_vmax = self._trf.transform([self.vmin, self.vmax])
   1710 if not np.isfinite([t_vmin, t_vmax]).all():
-> 1711     raise ValueError("Invalid vmin or vmax")
   1712 t_value -= t_vmin
   1713 t_value /= (t_vmax - t_vmin)

ValueError: Invalid vmin or vmax

<Figure size 2000x400 with 6 Axes>

In [51]:

stats_and_plot([2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019], [2015, 2016, 2017], 
               [], [2011, 2012, 2013, 2014, 2015, 2016], [], [],
               idepth, jdepth, Haro_Boundary, 'Haro_Boundary', 'alldepths')

In [52]:

stats_and_plot([2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019], [2015, 2016, 2017], 
               [], [2011, 2012, 2013, 2014, 2015, 2016], [], [],
               idepth, jdepth, JdF_west, 'JdF_west', 'alldepths')

In [53]:

stats_and_plot([2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019], [2015, 2016, 2017], 
               [], [2011, 2012, 2013, 2014, 2015, 2016], [], [],
               idepth, jdepth, JdF_east, 'JdF_east', 'alldepths')

Above 5 m¶

In [ ]:

idepth = 0
jdepth = 5

number, results = stats_allregions([2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019], 
                           [], [2015, 2016, 2017], [2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019], [], 
                           [],
                           idepth, jdepth, 'top5m')
for boxname in boxnames:
    if number > 0:
        print(f'\033[1m{boxname}\033[1m')
        display(results[results.region == boxname].loc[:, results.columns != 'region'].style.format(myformat).apply(highlight_max_min))

SS_all

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	6685	28.8	75.9	70.2	0.347	0.42	1.23e+02	42.3	53.9	0.405
202111	1644	-6.78	43.5	42.9	0.182	0.311	78.1	20.1	33.9	0.4

SoG_south

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	100	-9.0	33.0	31.7	0.126	-0.0767	-27.9	9.7	32.8	0.0501
202111	100	3.5	26.8	26.6	0.11	0.262	68.2	12.9	21.1	0.498

SoG_center

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	204	-13.8	48.7	46.7	0.335	0.482	1.31e+02	29.7	32.7	0.532
202111	204	-8.21	46.6	45.8	0.363	0.644	1.85e+02	36.0	24.9	0.708

SoG_north

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	145	-6.33	62.5	62.2	0.342	-0.113	-35.5	9.12	65.6	-0.033
202111	145	-1.51	56.2	56.2	0.343	0.311	78.5	18.1	46.8	0.33

PS_main

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	3285	45.2	89.6	77.4	0.431	0.498	1.5e+02	58.3	57.5	0.348

Haro_Boundary

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	99	2.16	36.4	36.3	0.188	-0.00352	1.4	2.16	36.4	-0.044
202111	99	6.91	34.3	33.6	0.208	0.314	74.4	15.9	26.8	0.384

JdF_west

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	414	42.1	78.8	66.6	0.387	0.523	1.48e+02	55.7	47.0	0.49
202111	130	2.66	26.0	25.9	0.127	0.116	31.0	5.32	24.1	0.156

JdF_east

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	617	10.4	36.0	34.4	0.157	0.0397	18.3	10.5	33.7	0.0668
202111	131	-3.88	20.3	19.9	0.0687	0.0868	15.5	5.27	18.8	0.201

Below 5 m¶

In [ ]:

idepth = 5
jdepth = 39

number, results = stats_allregions([2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019], 
                           [2015, 2016, 2017], [], [2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019], [], 
                           [],
                           idepth, jdepth, 'below5m')
for boxname in boxnames:
    if number > 0:
        print(f'\033[1m{boxname}\033[1m')
        display(results[results.region == boxname].loc[:, results.columns != 'region'].style.format(myformat).apply(highlight_max_min))

SS_all

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	9116	-21.8	41.0	34.7	0.113	-0.116	-43.0	22.8	36.2	-0.222
202111	9116	-2.58	29.4	29.2	0.0703	0.086	13.1	5.68	27.6	0.138

SoG_south

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	627	-12.3	26.4	23.3	0.081	-0.305	-72.8	17.6	23.4	-0.697
202111	627	1.52	16.1	16.0	0.0382	-0.024	-3.26	1.81	16.2	-0.0999

SoG_center

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	1520	-34.7	47.2	32.0	0.152	-0.218	-72.9	36.6	33.5	-0.582
202111	1520	0.29	25.1	25.1	0.0602	0.0691	12.4	3.74	24.0	0.0258

SoG_north

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	832	-34.1	50.9	37.7	0.211	-0.0065	-35.3	34.1	37.9	-0.054
202111	832	0.18	36.8	36.8	0.161	0.314	56.3	17.0	28.6	0.415

PS_main

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread

Haro_Boundary

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	615	-3.14	28.2	28.0	0.0928	-0.0395	-10.5	3.65	28.5	-0.0543
202111	615	1.13	27.0	26.9	0.0916	0.0608	12.5	3.09	26.0	0.105

JdF_west

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	816	-2.15	24.7	24.6	0.0342	-0.0803	-14.6	5.65	25.0	-0.112
202111	816	-3.07	21.4	21.2	0.0261	-0.0568	-11.9	4.81	21.5	-0.0763

JdF_east

	number	bias	rmse	crmse	swillmott	slopedev	const	systematic	nonsystematic	spread
201905R	625	-3.14	19.2	18.9	0.0348	-0.00167	-3.44	3.14	18.9	0.00956
202111	625	-2.52	17.3	17.1	0.0298	0.0507	6.65	3.65	16.5	0.112

In [ ]: