Figure 6

This notebook produces all subpanels of Figure 6, other auxillary figures, and the related statistical analysis.

import seaborn as sns
import pandas as pd
import matplotlib
from scipy import stats

Get data

We load the data produced by the previous workflow stored as ../data/results/pej0_GrIS_classification.csv. See Fig5.ipynb for more details.

df = pd.read_csv('../data/results/pej0_GrIS_classification.csv')
df
$\frac{P_e}{\ell}$ (m$^{-1}$) $J_0$ (m yr$^{-1}$) Speed diff (m yr$^{-1}$) Speed diff (%) Advance (km)
0 0.000021 -915.378881 3029.186179 66.610459 -17.0
1 0.000306 -200.805053 511.123464 84.192696 -2.3
2 -0.000047 -82.857761 456.284207 55.293407 -3.0
3 -0.000034 -450.936824 -435.655180 -12.603562 -0.3
4 0.000078 -1243.208582 -349.572187 -8.164585 -0.2
... ... ... ... ... ...
99 0.000234 -112.716518 613.747370 34.516832 -0.9
100 0.000044 -167.859915 1023.168837 55.202543 -4.1
101 0.000054 -1081.388452 2774.927415 51.468309 -8.8
102 0.000315 5.479272 467.138448 31.283371 -0.9
103 0.000231 -222.535111 573.128198 41.697811 -4.1

104 rows × 5 columns

This histogram gives a quick idea about the distribution of speed change in both absolute difference (the first histogram) and percentage difference (the second histogram). We use the absolute difference for the following analysis because it shows many outliers away from the peaked distribution centered at around 0 m/yr.

df.hist(column='Speed diff (m yr$^{-1}$)', bins=50);
../_images/Fig6_5_0.png
df.hist(column='Speed diff (%)', bins=50);
../_images/Fig6_6_0.png

Figure 6A

Visualization

We add a new column called Legend and classify the glaciers based on a simple threshold of speed change (here set to ±300 m/yr):

newclass = []
count = 0
for idx, row in df.iterrows():
    if row['Speed diff (m yr$^{-1}$)'] >= 300 or row['Speed diff (m yr$^{-1}$)'] <= -300:
        newclass.append(r'$|\Delta V| \geq$ 300 m yr$^{-1}$')
        count += 1
    else:
        newclass.append(r'$|\Delta V| <$ 300 m yr$^{-1}$')

df['Legend'] = newclass
print('How many glaciers have a higher speed change?: {}'.format(count))
print('How many glaciers have a lower speed change?:  {}'.format(len(df) - count))
How many glaciers have a higher speed change?: 54
How many glaciers have a lower speed change?:  50

Now we plot Figure 6A:

matplotlib.rc('font', size=20)
matplotlib.rc('axes', linewidth=2)

cmap = {r'$|\Delta V| \geq$ 300 m yr$^{-1}$': '#C94337', r'$|\Delta V| <$ 300 m yr$^{-1}$': '#8EB7E4'}
f = sns.jointplot(data=df, x=r'$\frac{P_e}{\ell}$ (m$^{-1}$)', y=r'$J_0$ (m yr$^{-1}$)', hue="Legend", palette=cmap, 
                  joint_kws={"s": 70, }, marginal_kws={'linewidth': 2}, height=8)
f.ax_joint.set_xlim([-0.0003, 0.0010])
f.ax_joint.set_ylim([-2000, 700])
f.ax_joint.set_xticks([-0.0003, 0, 0.0003, 0.0006, 0.0009])
f.ax_joint.tick_params(width=2, length=5)
f.ax_marg_x.tick_params(width=2, length=5)
f.ax_marg_y.tick_params(width=2, length=5)
legend = f.ax_joint.legend(loc='lower right')
legend.get_frame().set_linewidth(2)
legend.get_frame().set_edgecolor('xkcd:gray')
f.savefig('../data/results/Fig6-1.pdf')
../_images/Fig6_10_0.png

Statistical analysis

We use the two-sample Kolmogorov-Smirnov test to determine if the two kernel distributions shown in Figure 6A are drawn from the same probability distribution. The test is one dimensional, so we need to perform this twice, one for each axis.

unstable_j0 = df[df['Legend'] == r'$|\Delta V| \geq$ 300 m yr$^{-1}$'][r'$J_0$ (m yr$^{-1}$)']
unstable_pe = df[df['Legend'] == r'$|\Delta V| \geq$ 300 m yr$^{-1}$'][r'$\frac{P_e}{\ell}$ (m$^{-1}$)']
stable_j0 = df[df['Legend'] == r'$|\Delta V| <$ 300 m yr$^{-1}$'][r'$J_0$ (m yr$^{-1}$)']
stable_pe = df[df['Legend'] == r'$|\Delta V| <$ 300 m yr$^{-1}$'][r'$\frac{P_e}{\ell}$ (m$^{-1}$)']

We use the kstest fuction from the scipy package and perform a two-sided Kolmogorov–Smirnov test. The p value is sufficiently small to reject the null hypothesis that both kernels are drawn from the same probability distribution under a 95% confidence level.

stats.kstest(unstable_j0, stable_j0)
KstestResult(statistic=0.3451851851851852, pvalue=0.0028360448053456055)
stats.kstest(unstable_pe, stable_pe)
KstestResult(statistic=0.3237037037037037, pvalue=0.006217578940975077)

Figure 6B

First, this scatter plot shows how speed change and terminal retreat correlate to each other, motivating the needs to isolate that forcing in Figures 6B-C:

df.plot(x='Speed diff (m yr$^{-1}$)', y='Advance (km)', kind='scatter')
<AxesSubplot:xlabel='Speed diff (m yr$^{-1}$)', ylabel='Advance (km)'>
../_images/Fig6_17_1.png

First, let’s select glaciers with retreat > 0.5 km and plot them in Figure 6B using the same way as Figure 6A.

df_1 = df.loc[df['Advance (km)'] < -0.5]

cmap = {r'$|\Delta V| \geq$ 300 m yr$^{-1}$': '#C94337', r'$|\Delta V| <$ 300 m yr$^{-1}$': '#8EB7E4'}
f = sns.jointplot(data=df_1, x=r'$\frac{P_e}{\ell}$ (m$^{-1}$)', y=r'$J_0$ (m yr$^{-1}$)', hue="Legend", palette=cmap, 
                  joint_kws={"s": 70, }, marginal_kws={'linewidth': 2}, height=8)
f.ax_joint.set_xlim([-0.0003, 0.0010])
f.ax_joint.set_ylim([-2000, 700])
f.ax_joint.set_xticks([-0.0003, 0, 0.0003, 0.0006, 0.0009])
f.ax_joint.tick_params(width=2, length=5)
f.ax_marg_x.tick_params(width=2, length=5)
f.ax_marg_y.tick_params(width=2, length=5)
legend = f.ax_joint.legend(loc='lower right')
legend.get_frame().set_linewidth(2)
legend.get_frame().set_edgecolor('xkcd:gray')
f.savefig('../data/results/Fig6-2.pdf')
../_images/Fig6_19_0.png

Below is the same statistical test for this subset. Separation of Pe/ is still significant, but J0 is not.

unstable_j0_1 = df_1[df_1['Legend'] == r'$|\Delta V| \geq$ 300 m yr$^{-1}$'][r'$J_0$ (m yr$^{-1}$)']
unstable_pe_1 = df_1[df_1['Legend'] == r'$|\Delta V| \geq$ 300 m yr$^{-1}$'][r'$\frac{P_e}{\ell}$ (m$^{-1}$)']
stable_j0_1 = df_1[df_1['Legend'] == r'$|\Delta V| <$ 300 m yr$^{-1}$'][r'$J_0$ (m yr$^{-1}$)']
stable_pe_1 = df_1[df_1['Legend'] == r'$|\Delta V| <$ 300 m yr$^{-1}$'][r'$\frac{P_e}{\ell}$ (m$^{-1}$)']
stats.kstest(unstable_j0_1, stable_j0_1)
KstestResult(statistic=0.24154589371980675, pvalue=0.3662958100582304)
stats.kstest(unstable_pe_1, stable_pe_1)
KstestResult(statistic=0.36231884057971014, pvalue=0.04984912204217817)

Figure 6C

Now let’s select glaciers with retreat < 0.5 km and plot them in Figure 6C using the same way as Figure 6A.

df_2 = df.loc[df['Advance (km)'] > -0.5]

cmap = {r'$|\Delta V| \geq$ 300 m yr$^{-1}$': '#C94337', r'$|\Delta V| <$ 300 m yr$^{-1}$': '#8EB7E4'}
f = sns.jointplot(data=df_2, x=r'$\frac{P_e}{\ell}$ (m$^{-1}$)', y=r'$J_0$ (m yr$^{-1}$)', hue="Legend", palette=cmap, 
                  joint_kws={"s": 70, }, marginal_kws={'linewidth': 2}, height=8)
f.ax_joint.set_xlim([-0.0003, 0.0010])
f.ax_joint.set_ylim([-2000, 700])
f.ax_joint.set_xticks([-0.0003, 0, 0.0003, 0.0006, 0.0009])
f.ax_joint.tick_params(width=2, length=5)
f.ax_marg_x.tick_params(width=2, length=5)
f.ax_marg_y.tick_params(width=2, length=5)
legend = f.ax_joint.legend(loc='lower right')
legend.get_frame().set_linewidth(2)
legend.get_frame().set_edgecolor('xkcd:gray')
f.savefig('../data/results/Fig6-3.pdf')
../_images/Fig6_25_0.png

Below is the same statistical test for this subset. Separation of J0 is now significant, but Pe/ is not.

unstable_j0_2 = df_2[df_2['Legend'] == r'$|\Delta V| \geq$ 300 m yr$^{-1}$'][r'$J_0$ (m yr$^{-1}$)']
unstable_pe_2 = df_2[df_2['Legend'] == r'$|\Delta V| \geq$ 300 m yr$^{-1}$'][r'$\frac{P_e}{\ell}$ (m$^{-1}$)']
stable_j0_2 = df_2[df_2['Legend'] == r'$|\Delta V| <$ 300 m yr$^{-1}$'][r'$J_0$ (m yr$^{-1}$)']
stable_pe_2 = df_2[df_2['Legend'] == r'$|\Delta V| <$ 300 m yr$^{-1}$'][r'$\frac{P_e}{\ell}$ (m$^{-1}$)']
stats.kstest(unstable_j0_2, stable_j0_2)
KstestResult(statistic=0.95, pvalue=1.8017409321613442e-05)
stats.kstest(unstable_pe_2, stable_pe_2)
KstestResult(statistic=0.4142857142857143, pvalue=0.25506683332770297)