Figure 7#
This notebook produces a combined KDE figure for both the abrupt and gradual regime transition examples associated with this folder.
import warnings
import pickle
from tqdm import tqdm
from collections import OrderedDict
import pyleoclim as pyleo
import numpy as np
import ammonyte as amt
import seaborn as sns
import matplotlib.pyplot as plt
import scipy as sp
from scipy.signal import find_peaks
from scipy.stats import gaussian_kde
#We suppress warnings for these notebooks for presentation purposes. Best practice is to not do this though.
import warnings
warnings.filterwarnings('ignore')
Load ODP record#
#Defining group lists for easy loading
group_names = ['ODP 925'] #,'ODP 927','ODP 929','ODP 846','ODP 849']
series_list = []
color_list = sns.color_palette('colorblind')
for name in group_names:
with open('./data/LR04cores_spec_corr/'+name[-3:]+'_LR04age.txt','rb') as handle:
lines = handle.readlines()
time = []
d18O = []
for x in lines:
line_time = float(format(float(x.decode().split()[1]),'10f'))
line_d18O = float(format(float(x.decode().split()[2]),'10f'))
#There is a discontinuity in 927 around 4000 ka, we'll just exclude it
if line_time <= 4000:
time.append(line_time)
d18O.append(line_d18O)
series = pyleo.Series(value=d18O,
time=time,
label=name,
time_name='Yr',
time_unit='ka',
value_name=r'$\delta^{18}O$',
value_unit='permil')
series_list.append(series)
max_time = min([max(series.time) for series in series_list])
min_time = max([min(series.time) for series in series_list])
ms = pyleo.MultipleSeries([series.slice((min_time,max_time)).interp() for series in series_list])
fig,ax = ms.stackplot(colors=color_list[:len(ms.series_list)],figsize=(8,2*len(ms.series_list)))
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[4], line 5
2 color_list = sns.color_palette('colorblind')
4 for name in group_names:
----> 5 with open('./data/LR04cores_spec_corr/'+name[-3:]+'_LR04age.txt','rb') as handle:
6 lines = handle.readlines()
7 time = []
File ~/miniconda3/envs/ammonyte/lib/python3.10/site-packages/IPython/core/interactiveshell.py:310, in _modified_open(file, *args, **kwargs)
303 if file in {0, 1, 2}:
304 raise ValueError(
305 f"IPython won't let you open fd={file} by default "
306 "as it is likely to crash IPython. If you know what you are doing, "
307 "you can use builtins' open."
308 )
--> 310 return io_open(file, *args, **kwargs)
FileNotFoundError: [Errno 2] No such file or directory: './data/LR04cores_spec_corr/925_LR04age.txt'
Load NGRIP record#
end_time=10000
NGRIP_lipd = pyleo.Lipd('./data/8k_ice/NGRIP.NGRIP.2004.lpd')
NGRIP_tso = NGRIP_lipd.to_LipdSeriesList()
NGRIP_series = NGRIP_tso[0].slice((0,end_time))
NGRIP_series.time_unit = 'ka'
Disclaimer: LiPD files may be updated and modified to adhere to standards
reading: NGRIP.NGRIP.2004.lpd
Finished read: 1 record
extracting paleoData...
extracting: NGRIP.NGRIP.2004
Created time series: 3 entries
Both age and year information are available, using age
Both age and year information are available, using age
Both age and year information are available, using age
/var/folders/5k/0y4jsz592qq0y78c_0ddgcpm0000gn/T/ipykernel_11923/2600103537.py:3: DeprecationWarning: The Lipd class is being deprecated and will be removed in Pyleoclim v1.0.0. Functionalities will instead be handled by the pyLipd package.
NGRIP_lipd = pyleo.Lipd('./data/8k_ice/NGRIP.NGRIP.2004.lpd')
/Users/alexjames/Documents/GitHub/Pyleoclim_util/pyleoclim/core/lipdseries.py:133: UserWarning: auto_time_params is not specified. Currently default behavior sets this to True. In a future release this will be changed to False.
super(LipdSeries, self).__init__(time=time, value=value, time_name=time_name,
NGRIP_series.plot()
(<Figure size 1000x400 with 1 Axes>,
<Axes: xlabel='Age [ka]', ylabel='d18O [permil]'>)
Defining functions#
# This function can be used to calculate when transitions occur. It isn't used in this notebook as the calculations were performed previously
def detect_transitions(series,transition_interval=None):
'''Function to detect transitions across a confidence interval
Parameters
----------
series : pyleo.Series, amt.Series
Series to detect transitions upon
transition_interval : list,tuple
Upper and lower bound for the transition interval
Returns
-------
transitions : list
Timing of the transitions of the series across its confidence interval
'''
series_fine = series.interp(step=1)
if transition_interval is None:
upper, lower = amt.utils.sampling.confidence_interval(series)
else:
upper, lower = transition_interval
above_thresh = np.where(series_fine.value > upper,1,0)
below_thresh = np.where(series_fine.value < lower,1,0)
transition_above = np.diff(above_thresh)
transition_below = np.diff(below_thresh)
upper_trans = series_fine.time[1:][np.diff(above_thresh) != 0]
lower_trans = series_fine.time[1:][np.diff(below_thresh) != 0]
full_trans = np.zeros(len(transition_above))
last_above = 0
last_below = 0
for i in range(len(transition_above)):
above = transition_above[i]
below = transition_below[i]
if above != 0:
if last_below+above == 0:
loc = int((i+below_pointer)/2)
full_trans[loc] = 1
last_below=0
last_above = above
above_pointer = i
if below != 0:
if last_above + below == 0:
loc = int((i+above_pointer)/2)
full_trans[loc] = 1
last_above=0
last_below = below
below_pointer = i
transitions = series_fine.time[1:][full_trans != 0]
return transitions
Plotting#
The calculations used to produce these results were all done separately, so we just load the pickle files with the results here. We used the same parameteres as those from Figures 3 and 5.
#Production figure
color_list = sns.color_palette('colorblind',4)
SMALL_SIZE = 34
MEDIUM_SIZE = 34
BIGGER_SIZE = 34
plt.rc('font', size=SMALL_SIZE) # controls default text sizes
plt.rc('axes', titlesize=MEDIUM_SIZE) # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
plt.rc('xtick', labelsize=MEDIUM_SIZE) # fontsize of the tick labels
plt.rc('ytick', labelsize=MEDIUM_SIZE) # fontsize of the tick labels
plt.rc('legend', fontsize=BIGGER_SIZE) # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
fig,ax = plt.subplots(figsize=(28,18),nrows=2)
fig.tight_layout(h_pad=5, w_pad=5)
axes = ax.ravel()
with open('./data/leloup_paillard_transitions.pkl','rb') as handle:
trans_res = pickle.load(handle)
noise_levels = [2,1.5,1,.5]
for idx,level in enumerate(noise_levels):
transition_hist = trans_res[level]
kde = gaussian_kde(transition_hist,bw_method=.06)
x_axis = np.linspace(transition_hist.min(),transition_hist.max(),10000)
evaluated = kde.evaluate(x_axis)
evaluated /= max(evaluated)
axes[0].plot(x_axis,evaluated,label=f'S/N: {level}',color=color_list[idx],linewidth=4)
axes[0].grid(False)
axes[0].set_xlabel('Age [kyr BP]')
axes[0].set_ylabel('Normalized Transition Density')
axes[0].axvline(1000,color='grey',linestyle='dashed',alpha=.5,linewidth=4,label='Transition Point')
axes[0].legend(bbox_to_anchor=(1.2,1),loc='upper right')
with open('./data/8k_transitions.pkl','rb') as handle:
trans_res = pickle.load(handle)
noise_levels = [4,3,2,1]
for idx,level in enumerate(noise_levels):
transition_hist = trans_res[level]
kde = gaussian_kde(transition_hist,bw_method=.06)
x_axis = np.linspace(transition_hist.min(),transition_hist.max(),10000)
evaluated = kde.evaluate(x_axis)
evaluated /= max(evaluated)
axes[1].plot(x_axis,evaluated,label=f'S/N: {level}',color=color_list[idx],linewidth=4)
# ax.axvspan(7200,7800,color='grey',alpha=.3,label='Positive Range')
# ax.axvspan(8800,9600,color='grey',alpha=.3)
axes[1].grid(False)
axes[1].set_xlabel('Age [yr BP]')
axes[1].set_ylabel('Normalized Transition Density')
axes[1].axvline(7800,color='grey',linestyle='dashed',alpha=.5,linewidth=4,label='Transition Point')
axes[1].axvline(8400,color='grey',linestyle='dashed',alpha=.5,linewidth=4)
axes[1].legend(bbox_to_anchor=(1.2,1),loc='upper right')
<matplotlib.legend.Legend at 0x2c134cdf0>
