Synthetic Data#
The purpose of this notebook is to create synthetic data for the purpose of comparison as used in the original publication and for sensitivity tests in the KDE plots.
Note: This notebook assumes the existence of pickle files that need to have been created previously. If you are running this notebook on your machine, make sure you’ve successfully run both of the notebooks in the Loading Data folder.
To create synthetic versions of our ensembles, we substitute the real data with AR(1) noise that has a synthetic event with a magnitude of our choosing.
import pickle
import os
import pyleoclim as pyleo
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib.transforms as transforms
from matplotlib.ticker import FormatStrFormatter
from pylipd.lipd import LiPD
with open('../../data/pickle/preprocessed_series_dict.pkl','rb') as handle:
preprocessed_series_dict = pickle.load(handle)
with open('../../data/pickle/preprocessed_ens_dict.pkl','rb') as handle:
preprocessed_ens_dict = pickle.load(handle)
lat_dict = {series.lat:series.label for series in preprocessed_series_dict.values()}
sort_index = np.sort(np.array(list(lat_dict.keys())))[::-1]
sort_label = [lat_dict[lat] for lat in sort_index]
preprocessed_series_dict = {label:preprocessed_series_dict[label] for label in sort_label} #Sort by latitude
preprocessed_ens_dict = {label:preprocessed_ens_dict[label] for label in sort_label} #Sort by latitude
# define the function for the spike
def spike(t,delta,A):
''' Function to create a spike
t : float
time
A : float
amplitude
delta : float
parameter for the curve'''
f=1/(2*(len(t)-1))
y = (A/np.arctan(1/delta))*np.arctan(np.sin(2*np.pi*t*f)/delta)
return y
# define a function to add a spike to a series
def add_spike(series,xstart,xend,A,method='smooth'):
'''Function to add a spike to a pyleoclim series
Parameters
----------
series : pyleo.Series;
the series to which the spike will be added
xstart : float
the starting year of the spike
xend : float
the ending year of the spike
A : float
the amplitude of the spike
method : str
"smooth": add a spike by spike;
otherwise, directly add values of A at each timestep'''
new_series = series.copy()
x = new_series.time
y = new_series.value
if method=='smooth':
y[(x>=xstart)&(x<=xend)] += spike(np.arange(0,sum((x>=xstart)&(x<=xend))),0.02,A)
else:
y[(x>=xstart)&(x<=xend)] += [(x>=xstart)&(x<=xend)]+np.full(xend-xstart+1,A)
new_series.value = y
return new_series
This step is quite similar to the one we took in the Load_Data_Pyleoclim notebook, though this time we’re going to be creating our own signal, rather than using real data. We will still use the real chronologies though.
data_path = '../../data/LiPD/full'
D = LiPD()
D.load_from_dir(data_path)
lipd_records = D.get_all_dataset_names()
synthetic_ens_dict = {}
synthetic_signal_dict = {}
sn_ratios = [.5,1,2]
holocene_bounds = (0,10000)
#Set a signal to noise ratio
for sn in sn_ratios:
synthetic_ens_dict[sn] = {}
synthetic_signal_dict[sn] = {}
for record in lipd_records:
D = LiPD()
D.load(os.path.join(data_path,f'{record}.lpd'))
df = D.get_timeseries_essentials().iloc[0]
series = pyleo.GeoSeries(
time = df['time_values'],
value=df['paleoData_values'],
time_name = 'Age',
time_unit = 'yrs BP',
value_name = r'$\delta^{18} O$',
value_unit = u'‰',
label=record,
lat = df['geo_meanLat'],
lon=df['geo_meanLon'],
archiveType='speleothem',
dropna=False,
verbose = False
)
processed_series = series.slice(holocene_bounds).interp().standardize().detrend(method='savitzky-golay')
#Fit AR1 model
g = pyleo.utils.ar1_fit(y=processed_series.value,t=processed_series.time)
#Generate surrogate values according to ar1 model and the unprocessed series time axis
surr_value = pyleo.utils.tsmodel.ar1_model(t=series.time,tau=g)
surr_series = pyleo.Series(series.time,surr_value,verbose=False)
surr_spike_series = add_spike(surr_series,3900,4100,np.std(surr_series.value)*sn)
synthetic_signal_dict[sn][record] = surr_spike_series.value
ens_df = D.get_ensemble_tables().iloc[0]
ens_series_list = []
processed_ens_series_list = []
for i in range(1000): #We know there are 1000 ensemble members
ens_series = pyleo.GeoSeries(
time = ens_df['ensembleVariableValues'].T[i],
value= surr_spike_series.value,
time_name = 'Age',
time_unit = 'yrs BP',
value_name = r'$\delta^{18} O$',
value_unit = u'‰',
label=record,
lat = df['geo_meanLat'],
lon=df['geo_meanLon'],
archiveType='speleothem',
dropna=False,
verbose=False
)
ens_series_list.append(ens_series)
synthetic_ens = pyleo.EnsembleSeries(ens_series_list)
synthetic_ens_dict[sn][record] = synthetic_ens
Loading 14 LiPD files
0%| | 0/14 [00:00<?, ?it/s]
14%|████████████ | 2/14 [00:00<00:02, 4.45it/s]
29%|████████████████████████ | 4/14 [00:00<00:01, 6.79it/s]
36%|██████████████████████████████ | 5/14 [00:00<00:01, 5.52it/s]
43%|████████████████████████████████████ | 6/14 [00:01<00:01, 5.49it/s]
50%|██████████████████████████████████████████ | 7/14 [00:01<00:01, 5.90it/s]
57%|████████████████████████████████████████████████ | 8/14 [00:01<00:01, 4.36it/s]
64%|██████████████████████████████████████████████████████ | 9/14 [00:01<00:01, 4.43it/s]
71%|███████████████████████████████████████████████████████████▎ | 10/14 [00:01<00:00, 4.84it/s]
79%|█████████████████████████████████████████████████████████████████▏ | 11/14 [00:02<00:00, 5.28it/s]
86%|███████████████████████████████████████████████████████████████████████▏ | 12/14 [00:02<00:00, 5.11it/s]
93%|█████████████████████████████████████████████████████████████████████████████ | 13/14 [00:02<00:00, 4.72it/s]
100%|███████████████████████████████████████████████████████████████████████████████████| 14/14 [00:02<00:00, 5.27it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 8.62it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 8.57it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 3.05it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 3.04it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 11.67it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 10.90it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.19it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.18it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 5.77it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 5.75it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 7.89it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 7.87it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 2.98it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 2.97it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.99it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.98it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 6.48it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 6.46it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 7.58it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 7.56it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 5.23it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 5.22it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.34it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.34it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 14.18it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 11.17it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 3.01it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 3.01it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 11.56it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 11.71it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.21it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.21it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 5.92it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 5.91it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 7.78it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 7.76it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 2.96it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 2.96it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.97it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.96it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 6.48it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 6.46it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 7.46it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 7.44it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 5.19it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 5.18it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.38it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.38it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 13.32it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 11.44it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 3.02it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 3.01it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 11.11it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 11.50it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.18it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.17it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 5.94it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 5.92it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 7.78it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 7.76it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 3.00it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 2.99it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.90it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.89it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 6.53it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 6.51it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 7.37it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 7.34it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 5.21it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 5.20it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.38it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 4.38it/s]
Loaded..
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 14.23it/s]
Loaded..
with open('../../data/pickle/synthetic_signal_dict.pkl','wb') as handle:
pickle.dump(synthetic_signal_dict,handle)
# Create a figure with a specified size
fig = plt.figure(figsize=(8, 16))
# Set up plot parameters
xlim = [0, 10000]
n_ts = len(preprocessed_ens_dict)
fill_between_alpha = 0.2
cmap = 'tab20'
labels = 'auto'
ylabel_fontsize = 12
spine_lw = 1.5
grid_lw = 0.5
label_x_loc = -0.15
v_shift_factor = 1
linewidth = 1.5
ax = {}
left = 0
width = 1
height = 1 / n_ts
bottom = 1
synthetic_ens = synthetic_ens_dict[2] #Set the signal to noise ratio to 2
synthetic_ens = {label:synthetic_ens[label] for label in sort_label} #Sort by latitude
# Create a color palette with the same number of colors as the length of synthetic_ens
colors = sns.color_palette('tab20', n_colors=len(synthetic_ens))
for idx, pair in enumerate(synthetic_ens.items()):
color = colors[idx]
label, ens = pair
bottom -= height * v_shift_factor
ax[idx] = fig.add_axes([left, bottom, width, height])
# Plot the ensemble envelope
ens.common_time(time_axis=preprocessed_series_dict[label].time, bounds_error=False).plot_envelope(ax=ax[idx], shade_clr=color, curve_clr=color)
# Set plot properties for the main axis
ax[idx].patch.set_alpha(0)
ax[idx].set_xlim(xlim)
time_label = 'Years (BP)'
value_label = '$\delta^{18} O$ [‰]'
ax[idx].set_ylabel(value_label, weight='bold', size=ylabel_fontsize)
# Add labels to the plot
trans = transforms.blended_transform_factory(ax[idx].transAxes, ax[idx].transData)
ax[idx].text(-.1, 0, label, horizontalalignment='right', transform=trans, color=color, weight='bold')
ylim = ax[idx].get_ylim()
ax[idx].set_yticks([ylim[0], 0, ylim[-1]])
ax[idx].yaxis.set_major_formatter(FormatStrFormatter('%.1f'))
ax[idx].grid(False)
# Set spine and tick properties based on index
if idx % 2 == 0:
ax[idx].spines['left'].set_visible(True)
ax[idx].spines['left'].set_linewidth(spine_lw)
ax[idx].spines['left'].set_color(color)
ax[idx].spines['right'].set_visible(False)
ax[idx].yaxis.set_label_position('left')
ax[idx].yaxis.tick_left()
else:
ax[idx].spines['left'].set_visible(False)
ax[idx].spines['right'].set_visible(True)
ax[idx].spines['right'].set_linewidth(spine_lw)
ax[idx].spines['right'].set_color(color)
ax[idx].yaxis.set_label_position('right')
ax[idx].yaxis.tick_right()
ylim_mag = max(ylim) - min(ylim)
offset = ylim_mag * .05
# Set additional plot properties
ax[idx].yaxis.label.set_color(color)
ax[idx].tick_params(axis='y', colors=color)
ax[idx].spines['top'].set_visible(False)
ax[idx].spines['bottom'].set_visible(False)
ax[idx].tick_params(axis='x', which='both', length=0)
ax[idx].set_xlabel('')
ax[idx].set_xticklabels([])
ax[idx].legend([])
xt = ax[idx].get_xticks()[1:-1]
for x in xt:
ax[idx].axvline(x=x, color='lightgray', linewidth=grid_lw, ls='-', zorder=-1)
ax[idx].axhline(y=0, color='lightgray', linewidth=grid_lw, ls='-', zorder=-1)
ax[idx].invert_xaxis()
ax[idx].invert_yaxis()
ax[idx].axvspan(4100, 3900, color='lightgrey', alpha=0.5)
# Set up the x-axis label at the bottom
bottom -= height * (1 - v_shift_factor)
ax[n_ts] = fig.add_axes([left, bottom, width, height])
ax[n_ts].set_xlabel(time_label)
ax[n_ts].spines['left'].set_visible(False)
ax[n_ts].spines['right'].set_visible(False)
ax[n_ts].spines['bottom'].set_visible(True)
ax[n_ts].spines['bottom'].set_linewidth(spine_lw)
ax[n_ts].set_yticks([])
ax[n_ts].patch.set_alpha(0)
ax[n_ts].set_xlim(xlim)
ax[n_ts].grid(False)
ax[n_ts].tick_params(axis='x', which='both', length=3.5)
xt = ax[n_ts].get_xticks()[1:-1]
for x in xt:
ax[n_ts].axvline(x=x, color='lightgray', linewidth=grid_lw, ls='-', zorder=-1)
ax[n_ts].invert_xaxis()
ax[n_ts].invert_yaxis()
