<jupyter_start><jupyter_text>Silicon above 15 m<jupyter_code>import datetime as dt
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import numpy as np
import pandas as pd
import statsmodels.api as sm
import xarray as xr
foramt = "{:.2}"
myformat = {'bias': foramt, 'rmse': foramt, 'swillmott': foramt, 'slopedev': foramt, 'const': foramt,
                         'systematic': foramt, 'nonsystematic':foramt, 
                    'spread': foramt}
with xr.open_dataset('/home/sallen/MEOPAR/grid/mesh_mask201702.nc') as mesh:
    deptht = mesh.gdept_1d[0].values
def bias(df, obs, mod):
    diffy = df[mod] - df[obs]
    return diffy.count(), diffy.mean()
def rmse(df, obs, mod):
    return (np.sqrt(((df[mod] - df[obs])**2).mean()))
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())
def slope_inter(df, obs, mod):
    X = df[obs]
    X2 = X**2
    y = df[mod]
    X = np.column_stack((X, X2))
    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 predictions, model.params['x2'], model.params['x1'], model.params['const'], systematic, nonsyst
def spread(df, obs, mod):
    return 1 - ((df[mod] - df[mod].mean())**2).mean() / ((df[obs] - df[obs].mean())**2).mean() 
def read_pieces(pieces):
    temp1 = pd.read_csv(pieces[0])
    for piece in pieces[1:]:
        nextpiece = pd.read_csv(piece)
        temp1 = pd.concat([temp1, nextpiece], ignore_index=True)
    return temp1
def plot_and_stats(temp1, name, idepth, jdepth):
    fig, ax = plt.subplots(1, 1, figsize=(6, 5))
    vmax = 90
    vmin = 0
    condition = (temp1.k >= idepth) & (temp1.k <= jdepth)
    print (deptht[idepth], deptht[jdepth])
    counts, xedges, yedges, color = ax.hist2d(temp1.Si[(temp1.k >= idepth) & (temp1.k <= jdepth)], 
                                              temp1.mod_silicon[(temp1.k >= idepth) & (temp1.k <= jdepth)], 
                                           bins=np.arange(vmin, vmax, 1), norm=LogNorm());
    fig.colorbar(color)
    number, tbias = bias(temp1[(temp1.k >= idepth) & (temp1.k <= jdepth)], 'Si', 'mod_silicon')
    trmse = rmse(temp1[(temp1.k >= idepth) & (temp1.k <= jdepth)], 'Si', 'mod_silicon')
    tswillmott = swillmott(temp1[temp1.k >=idepth], 'Si', 'mod_silicon')
    predictions, m2, m, c, syst, nonsyst = slope_inter(temp1[(temp1.k >= idepth) & (temp1.k <= jdepth)], 'Si', 'mod_silicon')
    tspread = spread(temp1[temp1.k >= idepth], 'Si', 'mod_silicon')
    ax.plot([vmin, vmax], [vmin, vmax], 'w-');
    xr = np.arange(vmin, vmax, 0.5)
    ax.plot(temp1.Si[condition], predictions, 'r.');
    sc = 5
    sh = 2*sc-1
    bot = 65
    top = bot + 2*sh 
    ax.arrow(5, bot, 0, sh-np.abs(tbias)/2, head_width=0.5*sc, head_length=0.2*sc, length_includes_head=True)
    ax.arrow(5, top, 0, -sh+np.abs(tbias)/2, head_width=0.5*sc, head_length=0.2*sc, length_includes_head=True)
    ax.arrow(9, bot, 0, sh-syst/2, head_width=0.5*sc, head_length=0.2*sc, length_includes_head=True)
    ax.arrow(9, top, 0, -sh+syst/2, head_width=0.5*sc, head_length=0.2*sc, length_includes_head=True)
    ax.arrow(13, bot, 0, sh-nonsyst/2, head_width=0.5*sc, head_length=0.2*sc, length_includes_head=True)
    ax.arrow(13, top, 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,
          'swillmott': tswillmott,
          'slopedev': 1-m,
          'const': c,
          'systematic': syst,
          'nonsystematic': nonsyst,
          'spread': tspread}
    ax.text(5-1.2, 48, 'bias', rotation=90)
    ax.text(9-1.2, 48, 'systematic', rotation=90)
    ax.text(13-1.2, 42, 'non-systematic', rotation=90)
    ax.set_title(f'{name}, Silicon between 0 and 15 m');
    dCp2 = pd.DataFrame(data=Cp2, index=[name])
    return dCp2

pieces = ('/home/sallen/202007/H201812/ObsModel_201812_Bio_20150101-20151231.csv',
         '/home/sallen/202007/H201812/ObsModel_201812_Bio_20160101-20161231.csv',
         '/home/sallen/202007/H201812/ObsModel_201812_Bio_20170101-20171231.csv',
         '/home/sallen/202007/H201812/ObsModel_201812_PSF_20150101-20151231.csv',
  #       '/home/sallen/202007/H201905/ObsModel_201812_PSF_20160101-20161231.csv',
  #       '/home/sallen/202007/H201905/ObsModel_201812_PSF_20170101-20171231.csv',
         '/home/sallen/202007/H201812/ObsModel_H201812_pug_20150101_20151231.csv',
         '/home/sallen/202007/H201812/ObsModel_H201812_pug_20160101_20161231.csv',
         '/home/sallen/202007/H201812/ObsModel_H201812_pug_20170101_20171231.csv',
         '/home/sallen/202007/H201812/ObsModel_H201812_hplc_20150101_20151231.csv',
         '/home/sallen/202007/H201812/ObsModel_H201812_hplc_20160101_20161231.csv',
         '/home/sallen/202007/H201812/ObsModel_H201812_hplc_20170101_20171231.csv')
temp1 = read_pieces(pieces)
temp1['Si'] = temp1.Si.fillna(value=temp1['Silicate [umol/L]'])
idepth = 0
jdepth = 14
d201812 = plot_and_stats(temp1, 'H201812', idepth, jdepth)
d201812.style.format(myformat)
pieces = ('/home/sallen/202007/H201905/ObsModel_201905_Bio_20150101-20151231.csv',
         '/home/sallen/202007/H201905/ObsModel_201905_Bio_20160101-20161231.csv',
         '/home/sallen/202007/H201905/ObsModel_201905_Bio_20170101-20171231.csv',
         '/home/sallen/202007/H201905/ObsModel_201905_PSF_20150101-20151231.csv',
  #       '/home/sallen/202007/H201905/ObsModel_201905_PSF_20160101-20161231.csv',
  #       '/home/sallen/202007/H201905/ObsModel_201905_PSF_20170101-20171231.csv',
         '/home/sallen/202007/H201905/ObsModel_H201905_pug_20150101_20151231.csv',
         '/home/sallen/202007/H201905/ObsModel_H201905_pug_20160101_20161231.csv',
         '/home/sallen/202007/H201905/ObsModel_H201905_pug_20170101_20171231.csv',
         '/home/sallen/202007/H201905/ObsModel_H201905_hplc_20150101_20151231.csv',
         '/home/sallen/202007/H201905/ObsModel_H201905_hplc_20160101_20161231.csv',
         '/home/sallen/202007/H201905/ObsModel_H201905_hplc_20170101_20171231.csv')
temp1 = read_pieces(pieces)
temp1['Si'] = temp1.Si.fillna(value=temp1['Silicate [umol/L]'])

d201905 = plot_and_stats(temp1, 'H201905', idepth, jdepth)
d201905.style.format(myformat)
pieces = ('/home/sallen/202007/202007C-p2/ObsModel_202007Cp2_Bio_20150101-20151231.csv',
         '/home/sallen/202007/202007C-p2/ObsModel_202007Cp2_Bio_20160101-20161231.csv',
         '/home/sallen/202007/202007C-p2/ObsModel_202007Cp2_Bio_20170101-20171231.csv',
         '/home/sallen/202007/202007C-p2/ObsModel_202007Cp2_PSF_20150101-20151231.csv',
 #        '/home/sallen/202007/202007C-p2/ObsModel_202007Cp2_PSF_20160101-20161231.csv',
 #        '/home/sallen/202007/202007C-p2/ObsModel_202007Cp2_PSF_20170101-20171231.csv',
         '/home/sallen/202007/202007C-p2/ObsModel_202007Cp2_PUG_20150101-20151231.csv',
         '/home/sallen/202007/202007C-p2/ObsModel_202007Cp2_PUG_20160101-20161231.csv',
         '/home/sallen/202007/202007C-p2/ObsModel_202007Cp2_PUG_20170101-20171231.csv',
         '/home/sallen/202007/202007C-p2/ObsModel_202007Cp2_hplc_20150101_20151231.csv',
         '/home/sallen/202007/202007C-p2/ObsModel_202007Cp2_hplc_20160101_20161231.csv',
         '/home/sallen/202007/202007C-p2/ObsModel_202007Cp2_hplc_20170101_20171231.csv')
temp1 = read_pieces(pieces)
temp1['Si'] = temp1.Si.fillna(value=temp1['Silicate [umol/L]'])
dCp2 = plot_and_stats(temp1, 'Cp2', idepth, jdepth)
dCp2.style.format(myformat)
pieces = ('/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_bot_20150101_20151231.csv',
          '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_bot_20160101_20161231.csv',
          '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_bot_20170101_20171231.csv',
         '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_psf_20150101_20150331.csv',
          '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_psf_20150401_20150630.csv',
          '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_PSF_20150701-20151231.csv',
#          '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_psf_20160101_20161231.csv',
#          '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_psf_20170101_20171231.csv',
         '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_pug_20150101_20150331.csv',
          '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_pug_20150401_20150630.csv',
          '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_PUG_20150701-20151231.csv',
          '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_pug_20160101_20161231.csv',
          '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_pug_20170101_20171231.csv',
         '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_hplc_20150101_20151231.csv',
         '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_hplc_20160101_20161231.csv',
         '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_hplc_20170101_20171231.csv')
temp1 = read_pieces(pieces)
temp1['Si'] = temp1.Si.fillna(value=temp1['Silicate [umol/L]'])
dCp3 = plot_and_stats(temp1, 'Cp3', idepth, jdepth)
dCp3.style.format(myformat)
pieces = ('/home/sallen/202007/202007D-again/ObsModel_202007D-again_Bio_20150101-20151231.csv',
         '/home/sallen/202007/202007D-again/ObsModel_202007D-again_Bio_20160101-20161231.csv',
         '/home/sallen/202007/202007D-again/ObsModel_202007D-again_Bio_20170101-20171231.csv',
         '/home/sallen/202007/202007D-again/ObsModel_202007D-again_PSF_20150101-20151231.csv',
  #       '/home/sallen/202007/202007D-again/ObsModel_202007D-again_PSF_20160101-20161231.csv',
 #         '/home/sallen/202007/202007D-again/ObsModel_202007D-again_PSF_20170101-20171231.csv',
          '/home/sallen/202007/202007D-again/ObsModel_202007D-again_pug_20150101_20151231.csv',
          '/home/sallen/202007/202007D-again/ObsModel_202007D-again_pug_20160101_20161231.csv',
          '/home/sallen/202007/202007D-again/ObsModel_202007D-again_pug_20170101_20171231.csv',
          '/home/sallen/202007/202007D-again/ObsModel_202007D-again_hplc_20150101_20151231.csv',
          '/home/sallen/202007/202007D-again/ObsModel_202007D-again_hplc_20160101_20161231.csv',
          '/home/sallen/202007/202007D-again/ObsModel_202007D-again_hplc_20170101_20171231.csv'
        )
temp1 = read_pieces(pieces)
temp1['Si'] = temp1.Si.fillna(value=temp1['Silicate [umol/L]'])
Dagain = plot_and_stats(temp1, 'Dagain', idepth, jdepth)
Dagain.style.format(myformat)
pieces = ('/home/sallen/202007/202007D-lowR/ObsModel_202007D-lowR_Bio_20150101-20151231.csv',
         '/home/sallen/202007/202007D-lowR/ObsModel_202007D-lowR_bot_20160101_20161231.csv',
         '/home/sallen/202007/202007D-lowR/ObsModel_202007D-lowR_bot_20170101_20171231.csv',
         '/home/sallen/202007/202007D-lowR/ObsModel_202007D-lowR_PSF_20150101-20151231.csv',
 #        '/home/sallen/202007/202007D-lowR/ObsModel_202007D-lowR_psf_20160101_20161231.csv',
 #        '/home/sallen/202007/202007D-lowR/ObsModel_202007D-lowR_psf_20170101_20171231.csv',
         '/home/sallen/202007/202007D-lowR/ObsModel_202007D-lowR_pug_20150101_20151231.csv',
         '/home/sallen/202007/202007D-lowR/ObsModel_202007D-lowR_pug_20160101_20161231.csv',
         '/home/sallen/202007/202007D-lowR/ObsModel_202007D-lowR_pug_20170101_20171231.csv',
         '/home/sallen/202007/202007D-lowR/ObsModel_202007D-lowR_hplc_20150101_20151231.csv',
         '/home/sallen/202007/202007D-lowR/ObsModel_202007D-lowR_hplc_20160101_20161231.csv',
         '/home/sallen/202007/202007D-lowR/ObsModel_202007D-lowR_hplc_20170101_20171231.csv')
temp1 = read_pieces(pieces)
temp1['Si'] = temp1.Si.fillna(value=temp1['Silicate [umol/L]'])
DlowR = plot_and_stats(temp1, 'D-lowR', idepth, jdepth)
DlowR.style.format(myformat)
pieces = ('/home/sallen/202007/202007F/ObsModel_202007F_bot_20150101_20151231.csv',
          '/home/sallen/202007/202007F/ObsModel_202007F_bot_20160101_20161231.csv',
          '/home/sallen/202007/202007F/ObsModel_202007F_bot_20170101_20171231.csv',
         '/home/sallen/202007/202007F/ObsModel_202007F_psf_20150101_20151231.csv',
#          '/home/sallen/202007/202007F/ObsModel_202007F_psf_20160101_20161231.csv',
#          '/home/sallen/202007/202007F/ObsModel_202007F_psf_20170101_20171231.csv',
         '/home/sallen/202007/202007F/ObsModel_202007F_pug_20150101_20151231.csv',
          '/home/sallen/202007/202007F/ObsModel_202007F_pug_20160101_20161231.csv',
          '/home/sallen/202007/202007F/ObsModel_202007F_pug_20170101_20171231.csv',
         '/home/sallen/202007/202007F/ObsModel_202007F_hplc_20150101_20151231.csv',
         '/home/sallen/202007/202007F/ObsModel_202007F_hplc_20160101_20161231.csv',
         '/home/sallen/202007/202007F/ObsModel_202007F_hplc_20170101_20171231.csv')
temp1 = read_pieces(pieces)
temp1['Si'] = temp1.Si.fillna(value=temp1['Silicate [umol/L]'])
modF = plot_and_stats(temp1, 'F', idepth, jdepth)
modF.style.format(myformat)
pieces = ('/home/sallen/202007/202007G-p1/ObsModel_202007Gp1_bot_20150101_20151231.csv',
          '/home/sallen/202007/202007G-p1/ObsModel_202007Gp1_bot_20160101_20161231.csv',
          '/home/sallen/202007/202007G-p1/ObsModel_202007Gp1_bot_20170101_20171231.csv',
         '/home/sallen/202007/202007G-p1/ObsModel_202007Gp1_psf_20150101_20151231.csv',
#bd          '/home/sallen/202007/202007G-p1/ObsModel_202007Gp1_psf_20160101_20161231.csv',
# bd         '/home/sallen/202007/202007G-p1/ObsModel_202007Gp1_psf_20170101_20171231.csv',
         '/home/sallen/202007/202007G-p1/ObsModel_202007Gp1_pug_20150101_20151231.csv',
          '/home/sallen/202007/202007G-p1/ObsModel_202007Gp1_pug_20160101_20161231.csv',
          '/home/sallen/202007/202007G-p1/ObsModel_202007Gp1_pug_20170101_20171231.csv',
         '/home/sallen/202007/202007G-p1/ObsModel_202007Gp1_hplc_20150101_20151231.csv',
         '/home/sallen/202007/202007G-p1/ObsModel_202007Gp1_hplc_20160101_20161231.csv',
         '/home/sallen/202007/202007G-p1/ObsModel_202007Gp1_hplc_20170101_20171231.csv'
         )
temp1 = read_pieces(pieces)
temp1['Si'] = temp1.Si.fillna(value=temp1['Silicate [umol/L]'])
modGp1 = plot_and_stats(temp1, 'Gp1', idepth, jdepth)
modGp1.style.format(myformat)
pieces = ('/home/sallen/202007/202007G-p2/ObsModel_202007Gp2_bot_20150101_20151231.csv',
          '/home/sallen/202007/202007G-p2/ObsModel_202007Gp2f0_bot_20160101_20161231.csv',
          '/home/sallen/202007/202007G-p2/ObsModel_202007Gp2f0_bot_20170101_20171231.csv',
         '/home/sallen/202007/202007G-p2/ObsModel_202007Gp2_psf_20150101_20151231.csv',
#          '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_psf_20160101_20161231.csv',
#          '/home/sallen/202007/202007C-p3/ObsModel_202007Cp3_psf_20170101_20171231.csv',
         '/home/sallen/202007/202007G-p2/ObsModel_202007Gp2_pug_20150101_20151231.csv',
          '/home/sallen/202007/202007G-p2/ObsModel_202007Gp2f0_pug_20160101_20161231.csv',
          '/home/sallen/202007/202007G-p2/ObsModel_202007Gp2f0_pug_20170101_20171231.csv',
         '/home/sallen/202007/202007G-p2/ObsModel_202007Gp2_hplc_20150101_20151231.csv',
         '/home/sallen/202007/202007G-p2/ObsModel_202007Gp2f0_hplc_20160101_20161231.csv',
         '/home/sallen/202007/202007G-p2/ObsModel_202007Gp2f0_hplc_20170101_20171231.csv')
temp1 = read_pieces(pieces)
temp1['Si'] = temp1.Si.fillna(value=temp1['Silicate [umol/L]'])
modGp2 = plot_and_stats(temp1, 'Gp2', idepth, jdepth)
modGp2.style.format(myformat)
pieces = ('/home/sallen/202007/202007H/ObsModel_202007H_bot_20150101_20151231.csv',
          '/home/sallen/202007/202007H/ObsModel_202007H_bot_20160101_20161231.csv',
          '/home/sallen/202007/202007H/ObsModel_202007H_bot_20170101_20171231.csv',
          '/home/sallen/202007/202007H/ObsModel_202007H_hplc_20150101_20151231.csv',
          '/home/sallen/202007/202007H/ObsModel_202007H_hplc_20160101_20161231.csv',
          '/home/sallen/202007/202007H/ObsModel_202007H_hplc_20170101_20171231.csv',
          '/home/sallen/202007/202007H/ObsModel_202007H_pug_20150101_20151231.csv',
          '/home/sallen/202007/202007H/ObsModel_202007H_pug_20160101_20161231.csv',
          '/home/sallen/202007/202007H/ObsModel_202007H_pug_20170101_20171231.csv',
            '/home/sallen/202007/202007H/ObsModel_202007H_psf_20150101_20151231.csv',
         )
temp1 = read_pieces(pieces)
temp1['Si'] = temp1.Si.fillna(value=temp1['Silicate [umol/L]'])
modH = plot_and_stats(temp1, 'H', idepth, jdepth)
modH.style.format(myformat)
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))]
alltogether = pd.concat([d201812, d201905, dCp2, dCp3, Dagain, DlowR, modF, modGp1, modGp2, modH], axis=0)
foramt = "{:.2}"
alltogether.style.format(myformat).apply(highlight_max_min)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14165465552032286, 'fraction_numerical': 0.26082444508503894, 'mean_word_length': 4.016194331983805, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '7479', 'license': 'Apache-2.0', 'n_tokens_mistral': 10298, 'n_tokens_neox': 7365, 'n_words': 638, 'path': 'notebooks/PhysTuning/SurfaceSilicon.ipynb', 'repo_name': 'SalishSeaCast/analysis-susan', 'size': '17345'}