Holocene Ice Window Increment#
import pyleoclim as pyleo
import matplotlib.pyplot as plt
import numpy as np
from tqdm import tqdm
import ammonyte as amt
import seaborn as sns
import pandas as pd
from pylipd.lipd import LiPD
/Users/alexjames/miniconda3/envs/ammonyte/lib/python3.10/site-packages/pandas/core/arrays/masked.py:60: UserWarning: Pandas requires version '1.3.6' or newer of 'bottleneck' (version '1.3.5' currently installed).
from pandas.core import (
#We suppress all warnings for these notebooks for presentation purposes. Best practice is to not do this though.
import warnings
warnings.filterwarnings('ignore')
color_list = sns.color_palette('colorblind')
lipd_path = '../data/8k_ice'
all_files = LiPD()
if __name__=='__main__':
all_files.load_from_dir(lipd_path)
record_names = all_files.get_all_dataset_names()
Loading 8 LiPD files
0%| | 0/8 [00:00<?, ?it/s]
75%|███████████████████████████████████████████████████████████████▊ | 6/8 [00:00<00:00, 55.46it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 8/8 [00:00<00:00, 59.14it/s]
Loaded..
series_list = []
# We specify the indices of interest in each dataframe by hand here
index_dict = {
'GRIP.GRIP.1992' : 'd18O',
'Renland.Johnsen.1992' : 'd18O',
'EDML.Stenni.2010' : 'bagd18O',
'EPICADomeC.Stenni.2010' : 'bagd18O',
'Vostok.Vimeux.2002' : 'temperature',
'GISP2.Grootes.1997' : 'd18O',
'NGRIP.NGRIP.2004' : 'd18O',
'TALDICE.Mezgec.2017' : 'd18O',
}
for record in record_names:
d = LiPD()
d.load(f'{lipd_path}/{record}.lpd')
df = d.get_timeseries_essentials()
row = df[df['paleoData_variableName']==index_dict[record]][df['time_variableName']=='age']
lat = row['geo_meanLat'].to_numpy()[0]
lon = row['geo_meanLon'].to_numpy()[0]
elevation = row['geo_meanElev'].to_numpy()[0]
value = row['paleoData_values'].to_numpy()[0]
value_name = row['paleoData_variableName'].to_numpy()[0]
value_unit = row['paleoData_units'].to_numpy()[0]
time = row['time_values'].to_numpy()[0]
time_unit = row['time_units'].to_numpy()[0]
time_name = row['time_variableName'].to_numpy()[0]
label = row['dataSetName'].to_numpy()[0]
geo_series = pyleo.GeoSeries(time=time,
value=value,
lat=lat,
lon=lon,
elevation=elevation,
time_unit=time_unit,
time_name=time_name,
value_name=value_name,
value_unit=value_unit,
label=label,
archiveType='ice')
series_list.append(geo_series)
geo_ms = pyleo.MultipleGeoSeries(series_list)
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 75.99it/s]
Loaded..
NaNs have been detected and dropped.
Time axis values sorted in ascending order
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 47.23it/s]
Loaded..
Time axis values sorted in ascending order
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 55.88it/s]
Loaded..
NaNs have been detected and dropped.
Time axis values sorted in ascending order
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 55.26it/s]
Loaded..
Time axis values sorted in ascending order
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 55.80it/s]
Loaded..
NaNs have been detected and dropped.
Time axis values sorted in ascending order
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 65.66it/s]
Loaded..
NaNs have been detected and dropped.
Time axis values sorted in ascending order
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 75.58it/s]
Loaded..
Time axis values sorted in ascending order
Loading 1 LiPD files
0%| | 0/1 [00:00<?, ?it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 65.28it/s]
Loaded..
NaNs have been detected and dropped.
Time axis values sorted in ascending order
ms_list = []
for series in geo_ms.series_list:
if series.lat > 0 or series.label == 'EPICADomeC.Stenni.2010':
series.time_unit = 'Years BP'
ms_list.append(series)
end_time = 10000
ms_ice = pyleo.MultipleSeries([series.slice((0,end_time)) for series in ms_list])
ms_ice.stackplot(colors=color_list[:len(ms_ice.series_list)])
(<Figure size 640x480 with 6 Axes>,
{0: <Axes: ylabel='d18O [permil]'>,
1: <Axes: ylabel='bagd18O [permil]'>,
2: <Axes: ylabel='d18O [permil]'>,
3: <Axes: ylabel='d18O [permil]'>,
4: <Axes: ylabel='d18O [permil]'>,
'x_axis': <Axes: xlabel='Age [Years BP]'>})
m = 12
w_incre_list = np.arange(2,9,2)
res_dict = {}
lp_dict = {}
for idx,series in enumerate(ms_ice.series_list):
name = series.label
print(f'Analyzing {name}')
for w_incre in tqdm(w_incre_list):
series_slice = series.slice((0,end_time))
td = amt.TimeEmbeddedSeries(series=series.slice((0,end_time)),m=m) #Tau is selected via first minimum of mutual info
eps = td.find_epsilon(eps=1,target_density=.05,tolerance=.01,verbose=False)
rm = eps['Output']
lp = rm.laplacian_eigenmaps(w_size=20,w_incre=w_incre)
if name not in res_dict:
res_dict[name] = []
if name not in lp_dict:
lp_dict[name] = []
lp_dict[name].append(lp)
Analyzing Renland.Johnsen.1992
0%| | 0/4 [00:00<?, ?it/s]
25%|█████████████████████▎ | 1/4 [00:00<00:02, 1.30it/s]
50%|██████████████████████████████████████████▌ | 2/4 [00:01<00:01, 1.68it/s]
75%|███████████████████████████████████████████████████████████████▊ | 3/4 [00:01<00:00, 2.01it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:01<00:00, 2.27it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:01<00:00, 2.02it/s]
Analyzing EPICADomeC.Stenni.2010
0%| | 0/4 [00:00<?, ?it/s]
25%|█████████████████████▎ | 1/4 [00:00<00:02, 1.21it/s]
50%|██████████████████████████████████████████▌ | 2/4 [00:01<00:01, 1.51it/s]
75%|███████████████████████████████████████████████████████████████▊ | 3/4 [00:01<00:00, 1.78it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:02<00:00, 2.02it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:02<00:00, 1.81it/s]
Analyzing GISP2.Grootes.1997
0%| | 0/4 [00:00<?, ?it/s]
25%|█████████████████████▎ | 1/4 [00:01<00:05, 1.85s/it]
50%|██████████████████████████████████████████▌ | 2/4 [00:03<00:03, 1.52s/it]
75%|███████████████████████████████████████████████████████████████▊ | 3/4 [00:04<00:01, 1.22s/it]
100%|█████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:04<00:00, 1.04s/it]
100%|█████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:04<00:00, 1.19s/it]
Analyzing GRIP.GRIP.1992
0%| | 0/4 [00:00<?, ?it/s]
25%|█████████████████████▎ | 1/4 [00:04<00:14, 4.71s/it]
50%|██████████████████████████████████████████▌ | 2/4 [00:08<00:08, 4.07s/it]
75%|███████████████████████████████████████████████████████████████▊ | 3/4 [00:11<00:03, 3.52s/it]
100%|█████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:13<00:00, 3.16s/it]
100%|█████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:13<00:00, 3.45s/it]
Analyzing NGRIP.NGRIP.2004
0%| | 0/4 [00:00<?, ?it/s]
25%|█████████████████████▎ | 1/4 [00:01<00:03, 1.18s/it]
50%|██████████████████████████████████████████▌ | 2/4 [00:01<00:01, 1.21it/s]
75%|███████████████████████████████████████████████████████████████▊ | 3/4 [00:02<00:00, 1.42it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:02<00:00, 1.68it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:02<00:00, 1.46it/s]
name = ms_list[0].label
fig,ax = plt.subplots(nrows=len(lp_dict[name]),figsize=(10,4*len(w_incre_list)),sharex=True)
fig.subplots_adjust(hspace=.5)
axes = ax.ravel()
for idx,series in tqdm(enumerate(lp_dict[name])):
series.confidence_fill_plot(ax=axes[idx],title=f'W_incre = {w_incre_list[idx]}')
axes[idx].set_ylabel('FI')
if idx < len(axes)-1:
axes[idx].set_xlabel('')
axes[idx].get_legend().remove()
fig.suptitle(name)
0it [00:00, ?it/s]
1it [00:02, 2.30s/it]
2it [00:04, 2.30s/it]
3it [00:06, 2.29s/it]
4it [00:09, 2.29s/it]
4it [00:09, 2.29s/it]
Text(0.5, 0.98, 'Renland.Johnsen.1992')
name = ms_list[1].label
fig,ax = plt.subplots(nrows=len(lp_dict[name]),figsize=(10,4*len(w_incre_list)),sharex=True)
fig.subplots_adjust(hspace=.5)
axes = ax.ravel()
for idx,series in tqdm(enumerate(lp_dict[name])):
series.confidence_fill_plot(ax=axes[idx],title=f'W_incre = {w_incre_list[idx]}')
axes[idx].set_ylabel('FI')
if idx < len(axes)-1:
axes[idx].set_xlabel('')
axes[idx].get_legend().remove()
fig.suptitle(name)
0it [00:00, ?it/s]
1it [00:02, 2.28s/it]
2it [00:04, 2.28s/it]
3it [00:06, 2.27s/it]
4it [00:09, 2.27s/it]
4it [00:09, 2.27s/it]
Text(0.5, 0.98, 'EPICADomeC.Stenni.2010')
name = ms_list[2].label
fig,ax = plt.subplots(nrows=len(lp_dict[name]),figsize=(10,4*len(w_incre_list)),sharex=True)
fig.subplots_adjust(hspace=.5)
axes = ax.ravel()
for idx,series in tqdm(enumerate(lp_dict[name])):
series.confidence_fill_plot(ax=axes[idx],title=f'W_incre = {w_incre_list[idx]}')
axes[idx].set_ylabel('FI')
if idx < len(axes)-1:
axes[idx].set_xlabel('')
axes[idx].get_legend().remove()
fig.suptitle(name)
0it [00:00, ?it/s]
1it [00:02, 2.35s/it]
2it [00:04, 2.34s/it]
3it [00:07, 2.34s/it]
4it [00:09, 2.33s/it]
4it [00:09, 2.33s/it]
Text(0.5, 0.98, 'GISP2.Grootes.1997')
name = ms_list[3].label
fig,ax = plt.subplots(nrows=len(lp_dict[name]),figsize=(10,4*len(w_incre_list)),sharex=True)
fig.subplots_adjust(hspace=.5)
axes = ax.ravel()
for idx,series in tqdm(enumerate(lp_dict[name])):
series.confidence_fill_plot(ax=axes[idx],title=f'W_incre = {w_incre_list[idx]}')
axes[idx].set_ylabel('FI')
if idx < len(axes)-1:
axes[idx].set_xlabel('')
axes[idx].get_legend().remove()
fig.suptitle(name)
0it [00:00, ?it/s]
1it [00:02, 2.55s/it]
2it [00:05, 2.53s/it]
3it [00:07, 2.52s/it]
4it [00:10, 2.53s/it]
4it [00:10, 2.53s/it]
Text(0.5, 0.98, 'GRIP.GRIP.1992')
name = ms_list[4].label
fig,ax = plt.subplots(nrows=len(lp_dict[name]),figsize=(10,4*len(w_incre_list)),sharex=True)
fig.subplots_adjust(hspace=.5)
axes = ax.ravel()
for idx,series in tqdm(enumerate(lp_dict[name])):
series.confidence_fill_plot(ax=axes[idx],title=f'W_incre = {w_incre_list[idx]}')
axes[idx].set_ylabel('FI')
if idx < len(axes)-1:
axes[idx].set_xlabel('')
axes[idx].get_legend().remove()
fig.suptitle(name)
0it [00:00, ?it/s]
1it [00:02, 2.18s/it]
2it [00:04, 2.16s/it]
3it [00:06, 2.17s/it]
4it [00:08, 2.15s/it]
4it [00:08, 2.16s/it]
Text(0.5, 0.98, 'NGRIP.NGRIP.2004')
