<jupyter_start><jupyter_text>Recalibration of SMP density coefficients*<NAME>, Environment and Climate Change Canada, 2020*The matched snow pit and SMP measurements from Part 1 are used to recalibrate the bilinear regression model of [Proksch et al. (2015)](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2014JF003266). We apply a K-Folds approach to minimize sampling bias while evaluating the OLS regression.<jupyter_code># Import community packages
import string
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
plt.rcParams["font.family"] = "Times New Roman"

import matplotlib.cm as cm
import pickle

from scipy import stats
from statsmodels.formula.api import ols
from sklearn.model_selection import KFold, StratifiedShuffleSplit

# Seed to replicate the paper result exactly
RANDOM_SEED = 2019

# Load comparison result from Part 1
result = pd.read_pickle('./output/smp_pit_filtered.pkl')
result['force_log'] = np.log(result['force_median'])<jupyter_output><empty_output><jupyter_text>K-Folds OLS methodOLS regression with 10 folds to minimize sampling bias.Model coefficients and skill are evaluated as the mean of all folds.<jupyter_code>k_fold = KFold(n_splits = 10, shuffle = True, random_state = RANDOM_SEED)
rmse = []; error = []; r = []; params = None

# Split the dataset into 10 roughly equal groups, 
# train on all but one test group
for train_idx, test_idx in k_fold.split(result):
    train = result.iloc[train_idx]
    test = result.iloc[test_idx]
    
    model_rho = ols("RHO ~ force_log + force_log * l", train).fit()
    predict_rho = model_rho.predict(exog=dict(force_log=test['force_log'], l=test['l']))
    rmse = np.append(rmse, np.sqrt(np.mean((predict_rho - test['RHO'])**2)))
    r = np.append(r,np.corrcoef(predict_rho, test['RHO'])[1][0])
    error = np.append(error, predict_rho - test['RHO'])
    
    if params is None:
        params = model_rho.params.values
    else:
        params =  np.vstack((params, model_rho.params.values))
# K20a evaluation stats (from k-folds method)
# Metrics presented as mean of all fold permutations

k20a_rmse = rmse.mean()
k20a_bias = error.mean()
k20a_r2 = r.mean()**2

print('K20a recalibration evaluation')
print('N: %i' % len(result))
print('RMSE: %0.1f' % k20a_rmse)
print('RMSE percent: %0.2f' % np.round(k20a_rmse/result.RHO.mean(),2))
print('bias: %0.1f' % k20a_bias)
print('r^2: %0.2f' % k20a_r2)
model_k20a_coeff = [np.round(params[:,0].mean(),2), np.round(params[:,1].mean(),2),
              np.round(params[:,3].mean(),2), np.round(params[:,2].mean(),2)]
var_coeffs = [np.round(params[:,0].std(),2), np.round(params[:,1].std(),2),
              np.round(params[:,3].std(),2), np.round(params[:,2].std(),2)]

# Unbiased coeffs
print(model_k20a_coeff)

# Save coeffs
filename = './output/density_k20a_coeffs.pkl'
pickle.dump(model_k20a_coeff, open(filename, 'wb'))
# Apply the new coeffs to estimate density
k20a_rho = model_k20a_coeff[0]+(model_k20a_coeff[1]*result['force_log'])+ \
           (model_k20a_coeff[2]*result['force_log']*result['l'])+ \
           model_k20a_coeff[3]*result['l']<jupyter_output><empty_output><jupyter_text>K-folds OLS with outliers removedOutliers were defined as SMP/Cutter comparisons where error > than the 95th quantile in the K20a recalibration.  We justify this in the paper in the context of the matching procedure which cannot be assumed to be perfect.  Note that this removes a small number of comparisons(n=10) which are not isolated to any one profile.<jupyter_code>#Remove outliers
result_lim = result.copy()
result_lim['f_l'] = (result_lim['l'])*result_lim['force_log']
result_lim['abs_error'] = np.abs(k20a_rho - result_lim['RHO'])
q_95 = result_lim['abs_error'].quantile(0.95)
result_lim = result_lim[result_lim['abs_error'] < q_95]
n_removed  = len(result) - len(result_lim)

print('Error threshold: %i kg m^-3' % q_95)
print('Data points removed: %i' % n_removed)<jupyter_output><empty_output><jupyter_text>Figure 5 with caption Comparison of the SMP regression parameters and corresponding snow density observations. Parameters include log-transformed median force (ln(F ̃)), microstructure length scale (L) and an interaction term (f ̃L). Relationships are separated by ice surface environment.'<jupyter_code># Code to construct Figure 5
f, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True, figsize=(15,5))

axs = [ax1, ax2, ax3]
for n, ax in enumerate(axs):
    ax.text(0.02, 0.92, '('+string.ascii_lowercase[n]+')', transform=ax.transAxes, 
            size=16, weight='bold')

axis_value_size = 14
axis_label_szie = 14
point_size = 15

ax1.tick_params(axis='both', which='major', labelsize=axis_value_size)
ax2.tick_params(axis='both', which='major', labelsize=axis_value_size)
ax3.tick_params(axis='both', which='major', labelsize=axis_value_size)

kws_myi = dict(color = 'deepskyblue', s = point_size)
kws_fyi = dict(color = 'black', s = point_size)

ax1.scatter(result_lim[result_lim['ice_type'] == 'f']['force_log'], 
            result_lim[result_lim['ice_type'] == 'f']['RHO'], 
            label='FYI', **kws_fyi)
ax1.scatter(result_lim[result_lim['ice_type'] == 'm']['force_log'], 
            result_lim[result_lim['ice_type'] == 'm']['RHO'], 
            label='MYI', **kws_myi)


ax2.scatter(result_lim[result_lim['ice_type'] == 'f']['l'], 
            result_lim[result_lim['ice_type'] == 'f']['RHO'], 
            label='FYI', **kws_fyi)
ax2.scatter(result_lim[result_lim['ice_type'] == 'm']['l'], 
            result_lim[result_lim['ice_type'] == 'm']['RHO'], 
            label='MYI', **kws_myi)

ax3.scatter(result_lim[result_lim['ice_type'] == 'f']['f_l'], 
            result_lim[result_lim['ice_type'] == 'f']['RHO'], 
            label='FYI', **kws_fyi)
ax3.scatter(result_lim[result_lim['ice_type'] == 'm']['f_l'], 
            result_lim[result_lim['ice_type'] == 'm']['RHO'], 
            label='MYI', **kws_myi)

ax1.set_ylabel('Snow pit density [kg m$\mathregular{^{-3}}$]', fontsize = axis_label_szie)
ax1.set_xlabel(r'$ln(\tilde{F})$ [N]', fontsize = axis_label_szie)
ax2.set_xlabel(r'$L$ [mm]', fontsize = axis_label_szie)
ax3.set_xlabel(r'$ln(\tilde{F}) \cdot L$ [N mm]', fontsize = axis_label_szie)

ax1.set_xlim(-4, 4)
ax2.set_xlim(0, 1.5)
ax3.set_xlim(-3, 3)

ax3.legend(fontsize = 12)
ax3.legend(handletextpad=0, fontsize=12)


f.savefig('./output/figures/Fig05_RegressionTerms_lowres.png', format='png')
f.savefig('./output/figures/Fig05_RegressionTerms_production.pdf', format='pdf', dpi = 300)
# Correlation between snow pit observed density and median force.L for all comparisons
print(np.round(np.corrcoef(result_lim.RHO, result_lim.force_log)[0][1],2))
# Correlation between snow pit observed density and median force, microstructural length scale and the interaction term
result_lim.groupby(['ice_type'])[['RHO','force_log','l','f_l']].corr()
# Correlation between snow pit observed density and median force, microstructural length scale and the interaction term
result_lim.groupby(['campaign','ice_type'])[['RHO','force_log','l','f_l']].corr()
k_fold = KFold(10, True, RANDOM_SEED)
rmse = []
error = []
r = []
params = None

for train_idx, test_idx in k_fold.split(result_lim):
    train = result_lim.iloc[train_idx]
    test = result_lim.iloc[test_idx]
    
    model_rho = ols("RHO ~ force_log + force_log * l + l", train).fit()
    predict_rho = model_rho.predict(exog=dict(force_log=test['force_log'], l=test['l']))
    rmse = np.append(rmse, np.sqrt(np.mean((predict_rho - test['RHO'])**2)))
    r = np.append(r,np.corrcoef(predict_rho, test['RHO'])[1][0])
    error = np.append(error, predict_rho - test['RHO'])
    
    if params is None:
        params = model_rho.params.values
    else:
        params =  np.vstack((params, model_rho.params.values))
# K20a evaluation stats (from kfolds method)
k20b_rmse = rmse.mean()
k20b_bias = error.mean()
k20b_r2 = r.mean()**2

print('K20b recalibration evaluation')
print('N: %i' % len(result_lim))
print('RMSE: %0.1f' % k20b_rmse)
print('RMSE percent: %0.2f' % np.round(k20b_rmse/result_lim.RHO.mean(),2))
print('bias: %0.1f' % k20b_bias)
print('r^2: %0.2f' % k20b_r2)
model_k20b_coeff = [np.round(params[:,0].mean(),2), np.round(params[:,1].mean(),2),
              np.round(params[:,3].mean(),2), np.round(params[:,2].mean(),2)]
var_coeffs = [np.round(params[:,0].std(),2), np.round(params[:,1].std(),2),
              np.round(params[:,3].std(),2), np.round(params[:,2].std(),2)]

# Unbiased coeffs
print(model_k20b_coeff)

# Save coeffs
filename = './output/density_k20b_coeffs.pkl'
pickle.dump(model_k20b_coeff, open(filename, 'wb'))
k20b_rho = model_k20b_coeff[0]+(model_k20b_coeff[1]*result_lim['force_log'])+ \
           (model_k20b_coeff[2]*result_lim['force_log']*result_lim['l'])+ \
           model_k20b_coeff[3]*result_lim['l']<jupyter_output><empty_output><jupyter_text>Figure 4 with caption Evaluation of the SMP density model parametrization of Proksch et al. (2015) (P15; Top) and recalibrated coefficients for snow on sea ice (K19b; Bottom). In both cases the model is evaluated against manual density cutter measurements of snow density.<jupyter_code># Code to construct Figure 4
# P15 statistics
p2015_rmse = np.sqrt(np.mean((result['mean_samp']-result['RHO'])**2))
p2015_bias = (result['error']).mean()
p2015_r2 = np.ma.corrcoef(result['mean_samp'],result['RHO'])[0, 1]**2
p2015_n = len(result['mean_samp'])
p2015_p = stats.pearsonr(result['mean_samp'],result['RHO'])[1]


# Fig constants
axis_value_size = 14
axis_label_size = 14
rho_bin_size = 20 #in kg m-3
line_start_a = 100
line_end_a = 700
line_end_b = 550
point_size = 8
common_bin = np.arange(line_start_a,line_end_a, rho_bin_size)


f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(12,10))
f.subplots_adjust(hspace=.25)

axs = [ax1, ax2, ax3, ax4]
for n, ax in enumerate(axs):
    ax.text(0.02, 0.9, '('+string.ascii_lowercase[n]+')', transform=ax.transAxes, 
            size=20, weight='bold')

ax1.set_xlim(line_start_a,line_end_a)
ax1.set_ylim(line_start_a,line_end_a)
ax2.set_xlim(line_start_a,line_end_a)
ax3.set_xlim(line_start_a,line_end_b)
ax3.set_ylim(line_start_a,line_end_b)
ax4.set_xlim(line_start_a,line_end_b)

ax3.set_yticks(np.arange(line_start_a, line_end_b+1, 100.0))
ax1.tick_params(axis='both', which='major', labelsize=axis_value_size)
ax2.tick_params(axis='both', which='major', labelsize=axis_value_size)
ax3.tick_params(axis='both', which='major', labelsize=axis_value_size)
ax4.tick_params(axis='both', which='major', labelsize=axis_value_size)


ax1.scatter(result['RHO'], result['mean_samp'], 
            s = point_size, color ='black', zorder = 1000)
ax1.plot([line_start_a, line_end_a], [line_start_a, line_end_a], 
         'k-', color = 'gray' ,alpha= 0.8, zorder = 500)

hist_kws = dict(density=True, bins=common_bin, histtype="stepfilled", linewidth=1.25)
ax2.ticklabel_format(axis='y',style='sci', scilimits=(1,5), useMathText=False)
ax2.hist(result['RHO'], alpha = 1, edgecolor="black", 
         color = 'grey', label = 'Pit', **hist_kws)
ax2.hist(result['mean_samp'], alpha = 0.6, edgecolor="black", 
         color = 'deepskyblue', label = 'SMP', **hist_kws)

ax1.set_ylabel('P15 density [kg m$\mathregular{^{-3}}$]', fontsize=axis_label_size)
ax1.set_xlabel('Snow pit density [kg m$\mathregular{^{-3}}$]', fontsize=axis_label_size)
ax2.set_ylabel('Probability density', fontsize=axis_label_size)
ax2.set_xlabel('Snow density [kg m$\mathregular{^{-3}}$]', fontsize=axis_label_size)
ax2.legend(edgecolor = 'black',  fontsize=axis_value_size)

ax3.scatter(result_lim['RHO'], k20b_rho, s = point_size, color ='black',  zorder = 1000)
ax3.plot([line_start_a, line_end_b], [line_start_a, line_end_b], 
         'k-', color = 'grey' ,alpha= 0.8,  zorder = 500)

hist_kws = dict(density=True, bins = common_bin, histtype= "stepfilled", linewidth = 1.25)
ax4.ticklabel_format(axis='y',style='sci', scilimits=(1,5), useMathText=False)

ax4.hist(result_lim['RHO'], alpha = 1, edgecolor="black", color = 'grey', label = 'Pit', **hist_kws)
ax4.hist(k20b_rho, alpha = 0.6, edgecolor="black", color = 'deepskyblue', label = 'SMP', **hist_kws)

ax3.set_ylabel('K20b density [kg m$\mathregular{^{-3}}$]',fontsize=axis_label_size)
ax3.set_xlabel('Snow pit density [kg m$\mathregular{^{-3}}$]',fontsize=axis_label_size)
ax4.set_ylabel('Probability density',fontsize=axis_label_size)
ax4.set_xlabel('Snow density [kg m$\mathregular{^{-3}}$]',fontsize=axis_label_size)

# Display stats
ax1.text(550, 150,'N: %i \nRMSE: %i \nR$^2$: %0.2f'%(p2015_n, p2015_rmse, p2015_r2), fontsize=12)
ax3.text(440, 130,'N: %i \nRMSE: %i \nR$^2$: %0.2f'%(len(result_lim), k20b_rmse, k20b_r2),  fontsize=12)

f.savefig('./output/figures/Fig04_ModelEval_lowres.png', format='png')
f.savefig('./output/figures/Fig04_ModelEval_production.pdf', format='pdf', dpi = 300)
# Error metrics
def rmse(data):
    return np.sqrt(np.mean(data**2))

result_lim['model_rho'] = k20b_rho
result_lim['abs_error'] = np.abs(k20b_rho-result_lim['RHO']).values
result_lim['error'] = (k20b_rho-result_lim['RHO']).values

# Error by layer type
np.round(result_lim.groupby('TYPE')['error'].apply(rmse)/result_lim.groupby('TYPE')['model_rho'].mean(),3)
np.round(result_lim.groupby('campaign')['error'].apply(rmse)/result_lim.groupby('campaign')['model_rho'].mean(),3)
# Overall error in %
np.round((rmse(result_lim['error'])/result_lim['model_rho'].mean()),3)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16959761549925484, 'fraction_numerical': 0.03636363636363636, 'mean_word_length': 5.700449326010983, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 23, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '106370', 'license': 'MIT', 'n_tokens_mistral': 5391, 'n_tokens_neox': 5081, 'n_words': 973, 'path': 'Part_2_Calibration.ipynb', 'repo_name': 'kingjml/SMP-Sea-Ice', 'size': '13423'}