#!/usr/bin/env python
# coding: utf-8

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


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


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


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


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


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


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


# In[ ]: