MC-PCA on non-detrended speleothem data#
This Jupyter Notebook provides an analysis of speleothem oxygen isotope data using the Monte Carlo Principal Component Analysis (MC-PCA) method.
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.
import pickle
from tqdm import tqdm
import pyleoclim as pyleo
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import matplotlib as mpl
import seaborn as sns
from matplotlib import cm
from matplotlib.gridspec import GridSpec
with open('../../data/pickle/ens_dict.pkl','rb') as handle: #Using un preprocessed data
ens_dict = pickle.load(handle)
#Separating spels into two sets. Exclude Mawmluh for PCA
EA_set = [
'Jiuxian.China.2010',
'Xianglong.China.2018',
'Heshang.China.2008',
'Lianhua.China.2013',
'Dongge.China.2005',
# 'Guizhouxinv.China.2023', this gets excluded due to its short length
'Jiulong.China.2021',
'Liuli.China.2021',
]
IB_set = [
'Sahiya.India.2017',
'Hoq.Yemen.2013',
'Tangga.Indonesia.2018',
'Qunf.Oman.2023',
'LaVierge.Rodrigues.2018'
]
all_set = EA_set + IB_set
First we’ll analyze all of the speleothems together:
# Create multipleensemblegeoseries object, slicing off each timeseries at 3000 and 5000 years bp
megs = pyleo.MulEnsGeoSeries(list(ens_dict[label].interp().slice((3000,5000)) for label in all_set),label="All Speleothems")
#Applying mcpca to the multipleensemblegeoseries object
megs_pca = megs.mcpca(seed=42)
Iterating over simulations: 0%| | 0/1000 [00:00<?, ?it/s]
Iterating over simulations: 2%|█▎ | 25/1000 [00:00<00:03, 246.62it/s]
Iterating over simulations: 5%|██▋ | 52/1000 [00:00<00:03, 258.81it/s]
Iterating over simulations: 8%|████▏ | 80/1000 [00:00<00:03, 264.67it/s]
Iterating over simulations: 11%|█████▌ | 108/1000 [00:00<00:03, 267.81it/s]
Iterating over simulations: 14%|██████▉ | 135/1000 [00:00<00:03, 217.84it/s]
Iterating over simulations: 16%|████████▎ | 163/1000 [00:00<00:03, 235.33it/s]
Iterating over simulations: 19%|█████████▋ | 190/1000 [00:00<00:03, 243.59it/s]
Iterating over simulations: 22%|███████████▏ | 219/1000 [00:00<00:03, 255.35it/s]
Iterating over simulations: 25%|████████████▌ | 247/1000 [00:00<00:02, 262.49it/s]
Iterating over simulations: 28%|██████████████ | 276/1000 [00:01<00:02, 269.14it/s]
Iterating over simulations: 31%|███████████████▌ | 306/1000 [00:01<00:02, 277.34it/s]
Iterating over simulations: 34%|█████████████████▏ | 336/1000 [00:01<00:02, 282.68it/s]
Iterating over simulations: 37%|██████████████████▋ | 366/1000 [00:01<00:02, 285.86it/s]
Iterating over simulations: 40%|████████████████████▏ | 396/1000 [00:01<00:02, 288.59it/s]
Iterating over simulations: 43%|█████████████████████▋ | 426/1000 [00:01<00:01, 289.13it/s]
Iterating over simulations: 46%|███████████████████████▎ | 456/1000 [00:01<00:01, 290.45it/s]
Iterating over simulations: 49%|████████████████████████▊ | 486/1000 [00:01<00:01, 290.38it/s]
Iterating over simulations: 52%|██████████████████████████▎ | 516/1000 [00:01<00:01, 292.21it/s]
Iterating over simulations: 55%|███████████████████████████▊ | 546/1000 [00:02<00:01, 288.19it/s]
Iterating over simulations: 57%|█████████████████████████████▎ | 575/1000 [00:02<00:01, 285.32it/s]
Iterating over simulations: 60%|██████████████████████████████▊ | 604/1000 [00:02<00:01, 283.32it/s]
Iterating over simulations: 63%|████████████████████████████████▎ | 633/1000 [00:02<00:01, 281.19it/s]
Iterating over simulations: 66%|█████████████████████████████████▊ | 662/1000 [00:02<00:01, 282.81it/s]
Iterating over simulations: 69%|███████████████████████████████████▎ | 692/1000 [00:02<00:01, 285.18it/s]
Iterating over simulations: 72%|████████████████████████████████████▊ | 722/1000 [00:02<00:00, 288.19it/s]
Iterating over simulations: 75%|██████████████████████████████████████▎ | 752/1000 [00:02<00:00, 289.43it/s]
Iterating over simulations: 78%|███████████████████████████████████████▉ | 782/1000 [00:02<00:00, 289.67it/s]
Iterating over simulations: 81%|█████████████████████████████████████████▎ | 811/1000 [00:02<00:00, 288.41it/s]
Iterating over simulations: 84%|██████████████████████████████████████████▊ | 840/1000 [00:03<00:00, 286.20it/s]
Iterating over simulations: 87%|████████████████████████████████████████████▎ | 869/1000 [00:03<00:00, 286.66it/s]
Iterating over simulations: 90%|█████████████████████████████████████████████▊ | 899/1000 [00:03<00:00, 288.55it/s]
Iterating over simulations: 93%|███████████████████████████████████████████████▍ | 929/1000 [00:03<00:00, 290.23it/s]
Iterating over simulations: 96%|████████████████████████████████████████████████▉ | 959/1000 [00:03<00:00, 290.65it/s]
Iterating over simulations: 99%|██████████████████████████████████████████████████▍| 989/1000 [00:03<00:00, 291.67it/s]
Iterating over simulations: 100%|██████████████████████████████████████████████████| 1000/1000 [00:03<00:00, 278.96it/s]
#plotting
fig,ax = megs_pca.screeplot(figsize=(8,8))
fig,ax = megs_pca.modeplot(index=0,scatter_kwargs={'sizes':(20,300),'linewidth':.1},map_kwargs={'projection':'Mollweide','extent':[20,150,-25,45]})
ax['pc'].set_xlim([5000,3000])
Performing spectral analysis on individual series: 0%| | 0/1000 [00:00<?, ?it/s]
Performing spectral analysis on individual series: 9%|██▍ | 86/1000 [00:00<00:01, 851.75it/s]
Performing spectral analysis on individual series: 18%|████▉ | 175/1000 [00:00<00:00, 868.82it/s]
Performing spectral analysis on individual series: 26%|███████▎ | 262/1000 [00:00<00:00, 851.01it/s]
Performing spectral analysis on individual series: 35%|█████████▋ | 348/1000 [00:00<00:00, 844.00it/s]
Performing spectral analysis on individual series: 43%|████████████ | 433/1000 [00:00<00:00, 839.51it/s]
Performing spectral analysis on individual series: 52%|██████████████▍ | 517/1000 [00:00<00:00, 836.54it/s]
Performing spectral analysis on individual series: 60%|████████████████▊ | 601/1000 [00:00<00:00, 833.06it/s]
Performing spectral analysis on individual series: 69%|███████████████████▏ | 686/1000 [00:00<00:00, 835.40it/s]
Performing spectral analysis on individual series: 77%|█████████████████████▋ | 773/1000 [00:00<00:00, 845.93it/s]
Performing spectral analysis on individual series: 86%|████████████████████████▏ | 862/1000 [00:01<00:00, 857.87it/s]
Performing spectral analysis on individual series: 95%|██████████████████████████▌ | 948/1000 [00:01<00:00, 851.63it/s]
Performing spectral analysis on individual series: 100%|███████████████████████████| 1000/1000 [00:01<00:00, 847.17it/s]
(5000.0, 3000.0)
fig,ax = megs_pca.modeplot(index=1,scatter_kwargs={'sizes':(20,300),'linewidth':.1},map_kwargs={'projection':'Mollweide','extent':[20,150,-25,45]})
ax['pc'].set_xlim([5000,3000])
Performing spectral analysis on individual series: 0%| | 0/1000 [00:00<?, ?it/s]
Performing spectral analysis on individual series: 9%|██▍ | 86/1000 [00:00<00:01, 852.59it/s]
Performing spectral analysis on individual series: 17%|████▊ | 172/1000 [00:00<00:00, 848.51it/s]
Performing spectral analysis on individual series: 26%|███████▏ | 258/1000 [00:00<00:00, 853.07it/s]
Performing spectral analysis on individual series: 34%|█████████▋ | 345/1000 [00:00<00:00, 858.52it/s]
Performing spectral analysis on individual series: 43%|████████████ | 433/1000 [00:00<00:00, 864.70it/s]
Performing spectral analysis on individual series: 52%|██████████████▌ | 522/1000 [00:00<00:00, 870.15it/s]
Performing spectral analysis on individual series: 61%|█████████████████ | 610/1000 [00:00<00:00, 868.22it/s]
Performing spectral analysis on individual series: 70%|███████████████████▌ | 699/1000 [00:00<00:00, 872.55it/s]
Performing spectral analysis on individual series: 79%|██████████████████████ | 788/1000 [00:00<00:00, 877.91it/s]
Performing spectral analysis on individual series: 88%|████████████████████████▌ | 877/1000 [00:01<00:00, 879.19it/s]
Performing spectral analysis on individual series: 97%|███████████████████████████ | 966/1000 [00:01<00:00, 881.80it/s]
Performing spectral analysis on individual series: 100%|███████████████████████████| 1000/1000 [00:01<00:00, 870.51it/s]
(5000.0, 3000.0)
Now just the East Asian speleothems:
megs_ea = pyleo.MulEnsGeoSeries(list(ens_dict[label].interp().slice((3000,5000)) for label in EA_set),label='East Asian Speleothems')
megs_ea_pca = megs_ea.mcpca(seed=42)
Iterating over simulations: 0%| | 0/1000 [00:00<?, ?it/s]
Iterating over simulations: 4%|██▎ | 44/1000 [00:00<00:02, 434.77it/s]
Iterating over simulations: 9%|████▋ | 90/1000 [00:00<00:02, 444.36it/s]
Iterating over simulations: 14%|██████▉ | 136/1000 [00:00<00:01, 451.30it/s]
Iterating over simulations: 18%|█████████▎ | 182/1000 [00:00<00:01, 452.69it/s]
Iterating over simulations: 23%|███████████▋ | 228/1000 [00:00<00:01, 449.71it/s]
Iterating over simulations: 28%|██████████████ | 275/1000 [00:00<00:01, 453.79it/s]
Iterating over simulations: 32%|████████████████▍ | 323/1000 [00:00<00:01, 460.88it/s]
Iterating over simulations: 37%|██████████████████▉ | 371/1000 [00:00<00:01, 464.82it/s]
Iterating over simulations: 42%|█████████████████████▎ | 419/1000 [00:00<00:01, 466.05it/s]
Iterating over simulations: 47%|███████████████████████▊ | 467/1000 [00:01<00:01, 467.76it/s]
Iterating over simulations: 52%|██████████████████████████▎ | 515/1000 [00:01<00:01, 468.82it/s]
Iterating over simulations: 56%|████████████████████████████▋ | 563/1000 [00:01<00:00, 471.57it/s]
Iterating over simulations: 61%|███████████████████████████████▏ | 611/1000 [00:01<00:00, 473.93it/s]
Iterating over simulations: 66%|█████████████████████████████████▌ | 659/1000 [00:01<00:00, 473.89it/s]
Iterating over simulations: 71%|████████████████████████████████████ | 707/1000 [00:01<00:00, 472.96it/s]
Iterating over simulations: 76%|██████████████████████████████████████▌ | 755/1000 [00:01<00:00, 474.20it/s]
Iterating over simulations: 80%|████████████████████████████████████████▉ | 803/1000 [00:01<00:00, 472.35it/s]
Iterating over simulations: 85%|███████████████████████████████████████████▍ | 851/1000 [00:01<00:00, 473.20it/s]
Iterating over simulations: 90%|█████████████████████████████████████████████▊ | 899/1000 [00:01<00:00, 470.11it/s]
Iterating over simulations: 95%|████████████████████████████████████████████████▎ | 947/1000 [00:02<00:00, 469.13it/s]
Iterating over simulations: 99%|██████████████████████████████████████████████████▋| 994/1000 [00:02<00:00, 468.25it/s]
Iterating over simulations: 100%|██████████████████████████████████████████████████| 1000/1000 [00:02<00:00, 465.70it/s]
fig,ax = megs_ea_pca.screeplot()
fig,ax = megs_ea_pca.modeplot(index=0,scatter_kwargs={'sizes':(20,300),'linewidth':.1},map_kwargs={'projection':'Mollweide','extent':[80,150,10,55]})
ax['pc'].set_xlim([5000,3000])
Performing spectral analysis on individual series: 0%| | 0/1000 [00:00<?, ?it/s]
Performing spectral analysis on individual series: 9%|██▌ | 87/1000 [00:00<00:01, 863.30it/s]
Performing spectral analysis on individual series: 18%|████▉ | 177/1000 [00:00<00:00, 882.12it/s]
Performing spectral analysis on individual series: 27%|███████▍ | 266/1000 [00:00<00:00, 878.80it/s]
Performing spectral analysis on individual series: 36%|█████████▉ | 357/1000 [00:00<00:00, 886.97it/s]
Performing spectral analysis on individual series: 45%|████████████▌ | 447/1000 [00:00<00:00, 888.50it/s]
Performing spectral analysis on individual series: 54%|███████████████ | 536/1000 [00:00<00:00, 881.79it/s]
Performing spectral analysis on individual series: 62%|█████████████████▌ | 625/1000 [00:00<00:00, 882.13it/s]
Performing spectral analysis on individual series: 71%|███████████████████▉ | 714/1000 [00:00<00:00, 879.75it/s]
Performing spectral analysis on individual series: 80%|██████████████████████▌ | 805/1000 [00:00<00:00, 886.80it/s]
Performing spectral analysis on individual series: 89%|█████████████████████████ | 894/1000 [00:01<00:00, 885.94it/s]
Performing spectral analysis on individual series: 98%|███████████████████████████▌| 983/1000 [00:01<00:00, 884.53it/s]
Performing spectral analysis on individual series: 100%|███████████████████████████| 1000/1000 [00:01<00:00, 882.73it/s]
(5000.0, 3000.0)
fig,ax = megs_ea_pca.modeplot(index=1,scatter_kwargs={'sizes':(20,300),'linewidth':.1},map_kwargs={'projection':'Mollweide','extent':[80,150,10,55]})
ax['pc'].set_xlim([5000,3000])
Performing spectral analysis on individual series: 0%| | 0/1000 [00:00<?, ?it/s]
Performing spectral analysis on individual series: 9%|██▌ | 90/1000 [00:00<00:01, 891.70it/s]
Performing spectral analysis on individual series: 18%|█████ | 180/1000 [00:00<00:00, 868.55it/s]
Performing spectral analysis on individual series: 27%|███████▍ | 267/1000 [00:00<00:00, 867.83it/s]
Performing spectral analysis on individual series: 35%|█████████▉ | 354/1000 [00:00<00:00, 862.91it/s]
Performing spectral analysis on individual series: 44%|████████████▎ | 441/1000 [00:00<00:00, 864.73it/s]
Performing spectral analysis on individual series: 53%|██████████████▊ | 528/1000 [00:00<00:00, 866.14it/s]
Performing spectral analysis on individual series: 62%|█████████████████▎ | 618/1000 [00:00<00:00, 876.82it/s]
Performing spectral analysis on individual series: 71%|███████████████████▊ | 708/1000 [00:00<00:00, 880.46it/s]
Performing spectral analysis on individual series: 80%|██████████████████████▎ | 797/1000 [00:00<00:00, 883.33it/s]
Performing spectral analysis on individual series: 89%|████████████████████████▊ | 886/1000 [00:01<00:00, 879.60it/s]
Performing spectral analysis on individual series: 98%|███████████████████████████▎| 975/1000 [00:01<00:00, 882.62it/s]
Performing spectral analysis on individual series: 100%|███████████████████████████| 1000/1000 [00:01<00:00, 875.14it/s]
(5000.0, 3000.0)
Now just the speleothems from around the Indian Ocean Basin:
megs_ib = pyleo.MulEnsGeoSeries(list(ens_dict[label].interp().slice((3000,5000)) for label in IB_set),label='Indian Ocean Basin Speleothems')
megs_ib_pca = megs_ib.mcpca(seed=42)
Iterating over simulations: 0%| | 0/1000 [00:00<?, ?it/s]
Iterating over simulations: 6%|███ | 58/1000 [00:00<00:01, 570.33it/s]
Iterating over simulations: 12%|██████ | 118/1000 [00:00<00:01, 584.79it/s]
Iterating over simulations: 18%|█████████ | 177/1000 [00:00<00:01, 585.15it/s]
Iterating over simulations: 24%|████████████ | 236/1000 [00:00<00:01, 585.43it/s]
Iterating over simulations: 30%|███████████████ | 295/1000 [00:00<00:01, 585.73it/s]
Iterating over simulations: 36%|██████████████████▏ | 356/1000 [00:00<00:01, 590.86it/s]
Iterating over simulations: 42%|█████████████████████▏ | 416/1000 [00:00<00:00, 591.86it/s]
Iterating over simulations: 48%|████████████████████████▎ | 476/1000 [00:00<00:00, 593.83it/s]
Iterating over simulations: 54%|███████████████████████████▍ | 537/1000 [00:00<00:00, 595.86it/s]
Iterating over simulations: 60%|██████████████████████████████▍ | 598/1000 [00:01<00:00, 597.78it/s]
Iterating over simulations: 66%|█████████████████████████████████▌ | 658/1000 [00:01<00:00, 598.12it/s]
Iterating over simulations: 72%|████████████████████████████████████▋ | 719/1000 [00:01<00:00, 599.78it/s]
Iterating over simulations: 78%|███████████████████████████████████████▋ | 779/1000 [00:01<00:00, 595.52it/s]
Iterating over simulations: 84%|██████████████████████████████████████████▊ | 840/1000 [00:01<00:00, 598.14it/s]
Iterating over simulations: 90%|█████████████████████████████████████████████▉ | 900/1000 [00:01<00:00, 588.09it/s]
Iterating over simulations: 96%|████████████████████████████████████████████████▉ | 959/1000 [00:01<00:00, 587.00it/s]
Iterating over simulations: 100%|██████████████████████████████████████████████████| 1000/1000 [00:01<00:00, 591.00it/s]
fig,ax = megs_ib_pca.screeplot()
fig,ax = megs_ib_pca.modeplot(index=0,scatter_kwargs={'sizes':(20,300),'linewidth':.1},map_kwargs={'projection':'Mollweide','extent':[30,110,-30,35]})
ax['pc'].set_xlim([5000,3000])
Performing spectral analysis on individual series: 0%| | 0/1000 [00:00<?, ?it/s]
Performing spectral analysis on individual series: 8%|██▎ | 79/1000 [00:00<00:01, 788.79it/s]
Performing spectral analysis on individual series: 16%|████▍ | 160/1000 [00:00<00:01, 797.66it/s]
Performing spectral analysis on individual series: 24%|██████▋ | 240/1000 [00:00<00:00, 797.95it/s]
Performing spectral analysis on individual series: 32%|█████████ | 323/1000 [00:00<00:00, 808.53it/s]
Performing spectral analysis on individual series: 40%|███████████▎ | 404/1000 [00:00<00:00, 805.57it/s]
Performing spectral analysis on individual series: 48%|█████████████▌ | 485/1000 [00:00<00:00, 804.12it/s]
Performing spectral analysis on individual series: 57%|███████████████▊ | 566/1000 [00:00<00:00, 799.62it/s]
Performing spectral analysis on individual series: 65%|██████████████████ | 646/1000 [00:00<00:00, 797.88it/s]
Performing spectral analysis on individual series: 73%|████████████████████▍ | 729/1000 [00:00<00:00, 805.43it/s]
Performing spectral analysis on individual series: 81%|██████████████████████▋ | 812/1000 [00:01<00:00, 810.73it/s]
Performing spectral analysis on individual series: 89%|█████████████████████████ | 894/1000 [00:01<00:00, 807.06it/s]
Performing spectral analysis on individual series: 98%|███████████████████████████▎| 977/1000 [00:01<00:00, 811.99it/s]
Performing spectral analysis on individual series: 100%|███████████████████████████| 1000/1000 [00:01<00:00, 805.75it/s]
(5000.0, 3000.0)
fig,ax = megs_ib_pca.modeplot(index=1,scatter_kwargs={'sizes':(20,300),'linewidth':.1},map_kwargs={'projection':'Mollweide','extent':[30,110,-30,35]})
ax['pc'].set_xlim([5000,3000])
Performing spectral analysis on individual series: 0%| | 0/1000 [00:00<?, ?it/s]
Performing spectral analysis on individual series: 8%|██▏ | 77/1000 [00:00<00:01, 767.21it/s]
Performing spectral analysis on individual series: 16%|████▍ | 159/1000 [00:00<00:01, 794.44it/s]
Performing spectral analysis on individual series: 24%|██████▋ | 239/1000 [00:00<00:00, 785.59it/s]
Performing spectral analysis on individual series: 32%|████████▉ | 320/1000 [00:00<00:00, 792.71it/s]
Performing spectral analysis on individual series: 40%|███████████▏ | 400/1000 [00:00<00:00, 792.82it/s]
Performing spectral analysis on individual series: 48%|█████████████▍ | 480/1000 [00:00<00:00, 793.23it/s]
Performing spectral analysis on individual series: 56%|███████████████▋ | 560/1000 [00:00<00:00, 793.74it/s]
Performing spectral analysis on individual series: 64%|█████████████████▉ | 640/1000 [00:00<00:00, 789.48it/s]
Performing spectral analysis on individual series: 72%|████████████████████▏ | 720/1000 [00:00<00:00, 790.43it/s]
Performing spectral analysis on individual series: 80%|██████████████████████▍ | 801/1000 [00:01<00:00, 795.41it/s]
Performing spectral analysis on individual series: 88%|████████████████████████▋ | 881/1000 [00:01<00:00, 796.66it/s]
Performing spectral analysis on individual series: 96%|██████████████████████████▉ | 964/1000 [00:01<00:00, 804.85it/s]
Performing spectral analysis on individual series: 100%|███████████████████████████| 1000/1000 [00:01<00:00, 793.20it/s]
(5000.0, 3000.0)
