Orbital-scale variability over the last 1.5 Ma

%load_ext autoreload
%autoreload 2

import pyleoclim as pyleo
from pangaeapy.pandataset import PanDataSet

# load the dataset
dsb = PanDataSet('10.1594/PANGAEA.951382')
print(dsb.title)
display(dsb.data.head())

dsb.data.columns

dsb.data['Foram benth δ18O']

tsb = pyleo.GeoSeries(time=dsb.data['Age'], value=dsb.data['Foram benth δ18O'],
                      lat = dsb.data['Latitude'].iloc[0], lon = dsb.data['Longitude'].iloc[0],
                      elevation = dsb.data['Elevation'].iloc[0], archiveType = 'MarineSediment', 
                      time_name='Age', time_unit='ka BP', label = 'Benthic foram composite',
                      value_name='$\delta^{18}$O', value_unit=u'‰')

fig,ax = tsb.interp().dashboard(spectral_kwargs={'method':'mtm'}, plt_kwargs={'invert_yaxis':True})

fig, ax = tsb.plot(figsize=(10, 4),linewidth=0.5)
ax.invert_yaxis() # d18O is traditionally inverted
fig, ax = pyleo.add_GTS(fig, ax, ranks=['Period', 'Epoch'], location='above')

dsp = PanDataSet('10.1594/PANGAEA.951386')
print(dsp.title)
display(dsp.data.head())

tsp = pyleo.GeoSeries(time=dsp.data['Age'], value=dsp.data['G. bulloides δ18O'],
                      lat = dsp.data['Latitude'].iloc[0], lon = dsp.data['Longitude'].iloc[0],
                      elevation = dsb.data['Elevation'].iloc[0], archiveType = 'MarineSediment', 
                      time_name='Age', time_unit='ka BP', label = 'Planktonic record (G. bulloides)',
                      value_name='$\delta^{18}$O', value_unit=u'‰')
fig, ax = tsp.plot(figsize=(10, 4),linewidth=0.5,invert_yaxis=True, legend = False, ylabel = "G. Bulloides $\delta^{18}$O")
fig, ax = pyleo.add_GTS(fig, ax, ranks=['Period', 'Epoch'], location='above')

ms = tsp & tsb
fig, ax = ms.stackplot()
ax[0].invert_yaxis()
ax[1].invert_yaxis()

res, stats , sign =  pyleo.utils.tsbase.resolution(tsp.time)
import matplotlib.pyplot as plt
plt.hist(res)
stats

coh = tsp.wavelet_coherence(tsb,method='wwz')

coh_sig = coh.signif_test(number=100)

coh_sig.dashboard()

from climlab import constants as const
from climlab.solar.orbital import OrbitalTable
from climlab.solar.insolation import daily_insolation
import numpy as np

kyears = np.linspace(-1500, 0, 1001 ) # specify time interval and units
orb = OrbitalTable.interp(kyear=kyears) # subset of orbital parameters for specified time
days = np.linspace(0, const.days_per_year, 365)

Qsb = daily_insolation(dsb.data['Latitude'].iloc[0], days, orb) # generate insolation  at IODP339-U1385 latitude
Qsb_jja = np.mean(Qsb[:,151:243], axis=1)  # Julian days 152-243 are JJA

ts_ins = pyleo.Series(value = Qsb_jja, time = -kyears,
                      value_name = 'Insolation', time_name = 'Age',
                       value_unit = '$W/m^{2}$', time_unit='ky BP', label='Insolation', verbose=False )

# This code plots the planktonic d18O Series (tsp) and the insolation Series (ts_ins) on the same figure with two y-axes. The left y-axis shows d18O (inverted, as customary), and the right y-axis shows insolation. It converts both Pyleoclim Series to pandas Series for straightforward plotting with matplotlib, aligns x-limits to their common time coverage, and combines legends from both axes.
# Auto-generated by PaleoPal CodeGenerationAgent with contextual search

# Plot planktonic \u03b418O (tsp) and insolation (ts_ins) on dual y-axes, with insolation on the right
import matplotlib.pyplot as plt
import pandas as pd

# Convert Pyleoclim Series to pandas Series (index=time, values=value)
#tsp_pd = tsp.to_pandas(paleo_style=False)
#ts_ins_pd = ts_ins.to_pandas(paleo_style=False)

# Extract numpy arrays for plotting
#x1 = tsp_pd.index.values
#y1 = tsp_pd.values

#x2 = ts_ins_pd.index.values
#y2 = ts_ins_pd.values

#DK: Going through Pandas is technically correct but right now Matplotlib cannot plot large datetime values

x1 = tsp.time
y1 = tsp.value

x2 = ts_ins.time
y2 = ts_ins.value

# Create figure and primary axis
fig, ax = plt.subplots(figsize=(10, 4))

# Left y-axis: planktonic \u03b418O
line1, = ax.plot(x1, y1, color="tab:blue", linewidth=0.8, label="G Bulloides \u03b418O")
ax.set_ylabel("G. Bulloides \u03b418O (\u2030)", color="tab:blue")
ax.tick_params(axis="y", labelcolor="tab:blue")
ax.invert_yaxis()

# Right y-axis: insolation
ax2 = ax.twinx()
line2, = ax2.plot(x2, y2, color="tab:orange", linewidth=0.8, label="Insolation (JJA)")
ax2.set_ylabel("Insolation (W/m^2)", color="tab:orange")
ax2.tick_params(axis="y", labelcolor="tab:orange")

# X-axis label and common limits
ax.set_xlabel("Age (ka BP)")
try:
    xmin = max(float(x1.min()), float(x2.min()))
    xmax = min(float(x1.max()), float(x2.max()))
    if xmin < xmax:
        ax.set_xlim(xmin, xmax)
except Exception:
    # If limits cannot be computed, fall back to matplotlib defaults
    pass

# Combined legend
lines = [line1, line2]
labels = [l.get_label() for l in lines]
ax.legend(lines, labels, loc="upper right")

# Cosmetic tweaks
ax.grid(True, axis="x", linestyle=":", alpha=0.3)
#plt.tight_layout()
plt.show()

# Computes WWZ-based wavelet coherence between existing Series `tsp` (planktonic d18O) and `ts_ins` (JJA insolation), restricts both to their common time overlap, runs a Monte Carlo significance test, produces a coherence dashboard plot, and returns coherence objects plus phase statistics for precession/obliquity/eccentricity bands.
# Auto-generated by PaleoPal CodeGenerationAgent with contextual search

import numpy as np
import pyleoclim as pyleo

# Run wavelet coherence between planktonic \u03b418O (tsp) and JJA insolation (ts_ins)
# Prefer WWZ since tsp is unevenly spaced; ts_ins is evenly spaced but WTC can handle both.

# Ensure time axes are ascending (some paleo series may be descending)
try:
    tsp_sorted = tsp.sort(verbose=False, ascending=True, keep_log=False)
except Exception:
    tsp_sorted = tsp

try:
    ts_ins_sorted = ts_ins.sort(verbose=False, ascending=True, keep_log=False)
except Exception:
    ts_ins_sorted = ts_ins

# Restrict to common overlap to avoid edge effects and unnecessary computation
try:
    tmin = max(np.nanmin(tsp_sorted.time), np.nanmin(ts_ins_sorted.time))
    tmax = min(np.nanmax(tsp_sorted.time), np.nanmax(ts_ins_sorted.time))
    tsp_olap = tsp_sorted.slice((tmin, tmax))
    ts_ins_olap = ts_ins_sorted.slice((tmin, tmax))
except Exception:
    tsp_olap = tsp_sorted
    ts_ins_olap = ts_ins_sorted

# Coherence settings: orbital bands of interest ~20-200 kyr => frequencies 1/200 to 1/20 1/kyr
freq_kwargs = {"fmin": 1/200, "fmax": 1/20, "nf": 36}

# Compute wavelet coherence
coh_ins = tsp_olap.wavelet_coherence(
    ts_ins_olap,
    method="wwz",
    freq="log",
    freq_kwargs=freq_kwargs,
    verbose=True
)

# Significance testing (increase number for research-grade, e.g., 200+)
coh_ins_sig = coh_ins.signif_test(number=50)

# Dashboard plot (XWT + WTC + phase arrows + global coherence)
fig = coh_ins_sig.dashboard()

# Optional: phase statistics for specific scale bands (in kyr)
# Example: precession (19-23 kyr), obliquity (35-45 kyr), eccentricity (80-120 kyr)
phase_precession = coh_ins.phase_stats(scales=[19, 23])
phase_obliquity = coh_ins.phase_stats(scales=[35, 45])
phase_eccentricity = coh_ins.phase_stats(scales=[80, 120])

# Collect key outputs in a dict for downstream use
results = {
    "coherence": coh_ins,
    "coherence_signif": coh_ins_sig,
    "phase_precession": phase_precession,
    "phase_obliquity": phase_obliquity,
    "phase_eccentricity": phase_eccentricity,
    "overlap_timespan": (float(np.nanmin(tsp_olap.time)), float(np.nanmax(tsp_olap.time))) if hasattr(tsp_olap, "time") else None,
}

results