Temperature Scaling Behavior over the Holocene

# import the necessary packages
# Deal with downloading the LiPD files
import os
import requests
import zipfile
from pathlib import Path

#To mamipulate LiPD Files
from pylipd.lipd import LiPD

# Analysis
import pyleoclim as pyleo
import numpy as np
import pandas as pd

# Visualization
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib.colors import ListedColormap, BoundaryNorm
from matplotlib.cm import ScalarMappable

# To save
import pickle

# Link to LiPD files
temp12_url = 'https://lipdverse.org/Temp12k/current_version/Temp12k1_1_0.zip'

# Download and extract
folder_path = Path("../data/Temp12k")
zip_path = folder_path.with_suffix('.zip')

# Step 1: Check if folder exists, create if not:
folder_path.mkdir(parents=True, exist_ok=True)
# Step 2: Check if folder is empy
if not any(folder_path.iterdir()):
    print("Folder is empty. Downloading data...")

    # Step 3: Download zip file
    response = requests.get(temp12_url)
    with open(zip_path, 'wb') as f:
        f.write(response.content)
    print(f"Downloaded ZIP to {zip_path}")

    # Step 4: Extract zip contents into folder
    with zipfile.ZipFile(zip_path, 'r') as zip_ref:
        zip_ref.extractall(folder_path)
    print(f"Extracted contents to {folder_path}")

# This code snippet creates a LiPD object and loads data from the specified directory "../data/Temp12k". It then prints the names of all the datasets loaded into the LiPD object to verify that the data has been successfully imported.
# Auto-generated by PaleoPal CodeGenerationAgent with contextual search

# Import LiPD
#from pylipd.lipd import LiPD

# Create a LiPD object and load the data from the directory
lipd = LiPD()
lipd.load_from_dir("../data/Temp12k")

# Verify loaded datasets by printing their names
#print(lipd.get_all_dataset_names())

# This code snippet loads timeseries essentials from the Temp12k database using a LiPD object. It first ensures the data is loaded from the specified directory and then retrieves the timeseries essentials using the get_timeseries_essentials method, storing the result in a DataFrame for examination.
# Auto-generated by PaleoPal CodeGenerationAgent with contextual search

#from pylipd.lipd import LiPD

# Create a LiPD object if not already created
#lipd = LiPD()

# Load data from the directory if not already loaded
#lipd.load_from_dir('../data/Temp12k')

# Get the timeseries essentials
# This will return a pandas DataFrame or similar structure with timeseries information
#try:
df = lipd.get_timeseries_essentials()
    
    # Print the DataFrame to review the essential information
    #print(df_timeseries_essentials)
#except Exception as e:
    #print(f'An error occurred while getting timeseries essentials: {e}')

display(df.head())

df_filt = df[df['paleoData_variableName']=='temperature']
print(f"We have {len(df_filt)} timeseries")

display(df_filt.head())

# This code creates a Pyleoclim MultipleGeoSeries object from each row in a filtered DataFrame (df_filt) using essential timeseries information. It iterates through the DataFrame rows, constructing a GeoSeries object for each row, which is then appended to a list. Finally, a MultipleGeoSeries object is created from this list, allowing for further paleoclimate analysis.
# Auto-generated by PaleoPal CodeGenerationAgent with contextual search

# Import necessary library
#import pyleoclim as pyleo

# Create a list to hold GeoSeries objects
ts_list = []

# Iterate over each row in the DataFrame containing filtered timeseries\# Create a GeoSeries object for each row
for _, row in df_filt.iterrows():
    if row['time_values'] is not None and row['paleoData_values'] is not None: #there seems to be missing entries in the database (added by DK)
        if np.all(np.isnan(row['time_values'])) == False and np.all(np.isnan(row['paleoData_values'])) == False: #check that not all the values are false (added by DK)
            ts = pyleo.GeoSeries(
                time=row['time_values'],
                value=row['paleoData_values'],
                time_name=row['time_variableName'],
                value_name=row['paleoData_variableName'],
                time_unit=row['time_units'],
                value_unit=row['paleoData_units'],
                lat=row['geo_meanLat'],
                lon=row['geo_meanLon'],
                elevation=row.get('geo_meanElev'),  # Use get to handle missing values
                observationType=row['paleoData_proxy'],
                archiveType=row['archiveType'],
                verbose=False,
                label=row['dataSetName'] + '_' + row['paleoData_variableName']
            )
            ts_list.append(ts)
        else:
            pass
    else:
        pass

# Create a MultipleGeoSeries object from the list of GeoSeries objects
geoms = pyleo.MultipleGeoSeries(ts_list, label='Temp12k', time_unit='years BP')

# Display or process the MultipleGeoSeries object
# print(multiple_geoseries)  # Uncomment to see details of the object

print(f"We have {len(geoms.series_list)} timeseries left")

# This code snippet uses the MultipleGeoSeries object 'geoms' to create a map of the records using the 'map' method in Pyleoclim. The records are displayed in a Robinson projection with a specified figure size and color-coded by 'archiveType'.
# Auto-generated by PaleoPal CodeGenerationAgent with contextual search

# This code maps the records contained in the 'geoms' which is a MultipleGeoSeries object.
# Ensure that the pylipd and pyleoclim packages are installed and available for import.

# Import necessary package for plotting geographically distributed data
#import pyleoclim as pyleo

# 'geoms' is assumed to be a MultipleGeoSeries object already created from the DataFrame 'df_filt'

# Map the records
fig, ax = geoms.map(projection='Robinson', figsize=[12, 6], hue='archiveType')

# This code performs spectral analysis on a standardized MultipleGeoSeries object using the Weighted Wavelet Z-Transform (WWZ) method with a log frequency vector. The standardized data, 'geoms_standardized', is analyzed and the resulting Power Spectral Density (PSD) is plotted.
# Auto-generated by PaleoPal CodeGenerationAgent with contextual search

# Import necessary packages
#import pyleoclim as pyleo
#import matplotlib.pyplot as plt

# Standardize the 'geoms' MultipleGeoSeries object before spectral analysis
#geoms_standardized = geoms.standardize(keep_log=False) #keep log is not in MultipleSeries
geoms_standardized = geoms.standardize()

# Perform spectral analysis using the WWZ method with a log frequency vector

psd_wwz = geoms_standardized.spectral(method='wwz', freq='log')

# Plot the resultant Power Spectral Density (PSD) (might be too many to plot)
#fig, ax = psd_wwz.plot(title='PSD using WWZ method with log frequency vector')
#plt.show()

#with open("../data/temp12k_scaling.pkl", "wb") as f:   # 'wb' = write binary
    #pickle.dump(psd_wwz, f)

#If needed, re-open the pickle file. 
with open("../data/temp12k_scaling.pkl", "rb") as f:   # 'rb' = read binary
    psd_wwz = pickle.load(f)

label = []
lat = []
lon = []
archive = []
amp =[]
period = []
unit = []

#grab the information

for item in psd_wwz.psd_list:
    label.append(item.timeseries.label)
    lat.append(item.timeseries.lat)
    lon.append(item.timeseries.lon)
    archive.append(item.timeseries.archiveType)
    amp.append(item.amplitude)
    period.append(1/item.frequency)
    unit.append(item.period_unit)

# Put it in a DataFrame

df_psd = pd.DataFrame({'Name':label,
                       'Latitude':lat,
                       'Longitude':lon,
                       'Archive':archive,
                       'Amplitude':amp,
                       'Period':period,
                       'Units':unit})

display(df_psd.head())

# --- 0) Seaborn style
sns.set_theme(style="whitegrid", context="talk")  # 'talk' or 'poster' increases fonts
plt.rcParams.update({
    "font.size": 14,
    "axes.labelsize": 16,
    "axes.titlesize": 18,
    "xtick.labelsize": 12,
    "ytick.labelsize": 12,
    "legend.fontsize": 12,
})

# --- 1) Ensure arrays, explode to long-form (aligned lengths required)
# If arrays are lists/np arrays per row, this works in pandas >= 1.3
df_long = df_psd.copy()

# Coerce to list for safety
df_long["Amplitude"] = df_long["Amplitude"].apply(lambda x: np.asarray(x).tolist())
df_long["Period"]    = df_long["Period"].apply(lambda x: np.asarray(x).tolist())

# Add a series id (one per original row) so each becomes one line after explode
df_long["series_id"] = np.arange(len(df_long))

# Explode both columns (lists must have same length per row)
df_long = df_long.explode(["Period", "Amplitude"], ignore_index=True)

# Make numeric + drop NaNs / non-finite
df_long["Period"]    = pd.to_numeric(df_long["Period"], errors="coerce")
df_long["Amplitude"] = pd.to_numeric(df_long["Amplitude"], errors="coerce")
df_long = df_long.replace([np.inf, -np.inf], np.nan).dropna(subset=["Period", "Amplitude", "Latitude"])

# --- 2) Build absolute-latitude bins every 10°
abs_lat = np.abs(df_long["Latitude"])
# Bins: [0,10), [10,20), ..., [80,90]; include_lowest=True puts 0° into first bin
bins = np.arange(0, 100, 10)  # 0..90 by 10
labels = [f"{b:02d}–{b+10:02d}°" for b in bins[:-1]]  # 0–10°, 10–20°, ...
df_long["lat_bin"] = pd.cut(abs_lat, bins=bins, right=False, include_lowest=True, labels=labels)

# Keep only rows that fell into a defined bin
df_long = df_long.dropna(subset=["lat_bin"])

# --- 3) Custom palette per your mapping
# Start from RdBu_r, reverse, then splice BuGn_r(1:4) into positions [3:6]
cmap_base = sns.color_palette("RdBu_r", 9)[::-1]           # list of 9 RGB tuples
cmap_green = sns.color_palette("BuGn_r")                   # default ~10 colors
cmap_base[3:6] = cmap_green[1:4]                           # replace indices 3,4,5
listed_cmap = ListedColormap(cmap_base)

# Map bins to integer indices 0..8 in order of labels
bin_to_idx = {lab: i for i, lab in enumerate(labels)}
df_long["lat_bin_idx"] = df_long["lat_bin"].map(bin_to_idx)

# --- 4) Plot: one line per series, colored by lat_bin (discrete)
fig, ax = plt.subplots(figsize=(10, 6))

# We'll loop to guarantee exact color mapping & line-by-line control, while still using seaborn aesthetics
for (sid, latbin_idx), g in df_long.groupby(["series_id", "lat_bin_idx"], sort=False):
    ax.plot(
        g["Period"].values,
        g["Amplitude"].values,
        lw=1.2,
        alpha=0.9,
        color=listed_cmap(latbin_idx) if pd.notna(latbin_idx) else "0.6",
    )

# Log scale on x/y
ax.set_xscale("log")
ax.set_yscale("log")
ax.invert_xaxis()

# Reduce limits on y axis
ax.set_ylim(1e-6, 1e6)

# Redo the x ticks
period_ticks = [2, 5, 10, 20, 100, 1000, 10000, 100000, 1000000]
period_ticklabels = ['2', '5', '10', '20', '100', '1 k', '10 k', '100 k', '1 m']

ax.set_xticks(period_ticks)
ax.set_xticklabels(period_ticklabels)

# Labels (try to pull a single units label if unique)
y_units = None
if "Units" in df_long and df_long["Units"].nunique(dropna=True) == 1:
    y_units = df_long["Units"].dropna().iloc[0]

ax.set_xlabel(f"Period ({y_units})" if y_units else "Period")
ax.set_ylabel("PSD")

# --- 5) Discrete colorbar for the 10° abs-lat bins (top ticks + top label)
bounds = np.arange(-0.5, 9.5, 1)  # boundaries: -0.5, 0.5, ..., 8.5
norm = BoundaryNorm(boundaries=bounds, ncolors=listed_cmap.N)

sm = ScalarMappable(norm=norm, cmap=listed_cmap)
sm.set_array([])

# ticks exactly at bin edges → labels 0,10,...,90
tick_positions = np.arange(-0.5, 9.5, 1)
tick_labels    = [str(i * 10) for i in range(10)]  # 0..90

cbar = fig.colorbar(
    sm, ax=ax, boundaries=bounds, ticks=tick_positions, orientation='vertical'
)
cbar.set_ticklabels(tick_labels)            # numbers at the edges
cbar.ax.tick_params(labelsize=12)

# label
cbar.ax.set_ylabel("Latitude (°N/°S)", fontsize=16)

# --- 6) Style grid and axes lines
# Make grid dashed
ax.grid(True, which='major', linestyle='--', linewidth=0.8, color='gray', alpha=0.6)
ax.grid(True, which='minor', linestyle='--', linewidth=0.5, color='gray', alpha=0.4)

# Make the axis spines black and a bit thicker
for spine in ax.spines.values():
    spine.set_color('black')
    spine.set_linewidth(1.2)

# --- 7) Nice touches
plt.tight_layout()
plt.show()

# This code calculates the scaling exponent (beta) for a Power Spectral Density (PSD) object using the beta estimation method provided by the Pyleoclim library. The calculation is performed over three frequency ranges to yield beta_S, beta_M, and beta_L, which represent scaling exponents over short, medium, and long periods, respectively. The results are stored in the variables beta_S, beta_M, and beta_L, and are printed to the console.
# Auto-generated by PaleoPal CodeGenerationAgent with contextual search

# Calculate the scaling exponent (beta) for the psd_wwz object
# Import necessary libraries
#import pyleoclim as pyleo

# Define frequency bounds for beta calculations
frequency_ranges = [
    {'name': 'S', 'fmin': 1/100, 'fmax': 1/10},
    {'name': 'M', 'fmin': 1/1000, 'fmax': 1/100},
    {'name': 'L', 'fmin': 1/10000, 'fmax': 1/1000}
]

# Initialize dictionary to store beta values
beta_values = {}

# Calculate beta for each frequency range
for freq_range in frequency_ranges:
    # Calculate beta using the beta_est method
    beta_est_res = psd_wwz.beta_est(fmin=freq_range['fmin'], fmax=freq_range['fmax'])
    
    # Store result in the dictionary using the frequency range name as the key
    beta_values[f"beta_{freq_range['name']}"] = beta_est_res

# Retrieve the calculated beta values
beta_S = beta_values['beta_S'].beta_est_res['beta']
beta_M = beta_values['beta_M'].beta_est_res['beta']
beta_L = beta_values['beta_L'].beta_est_res['beta']

# Print the results
print(f"beta_S: {beta_S}")
print(f"beta_M: {beta_M}")
print(f"beta_L: {beta_L}")

# Transform into numpy arrays
beta_S = np.array(beta_S)
beta_M = np.array(beta_M)
beta_L = np.array(beta_L)

print(f"There are {len(beta_S[~np.isnan(beta_S)])} records for which the slope could be calculated in for periods between 10-100 year periods")
print(f"There are {len(beta_M[~np.isnan(beta_M)])} records for which the slope could be calculated in for periods between 100-1,000 year periods")
print(f"There are {len(beta_L[~np.isnan(beta_L)])} records for which the slope could be calculated in for periods between 1,000-10,000 year periods")

# --- Clean data: drop NaNs
beta_S_clean = np.array(beta_S)[~np.isnan(beta_S)]
beta_M_clean = np.array(beta_M)[~np.isnan(beta_M)]
beta_L_clean = np.array(beta_L)[~np.isnan(beta_L)]

# --- Compute medians and sample sizes
med_s = np.median(beta_S_clean)
med_m = np.median(beta_M_clean)
med_l = np.median(beta_L_clean)

n_s = beta_S_clean.shape[0]
n_m = beta_M_clean.shape[0]
n_l = beta_L_clean.shape[0]

# --- Plot
sns.set_theme(style="whitegrid", context="talk")
fig, ax = plt.subplots(figsize=[8, 6])

ax.set_xlim([-5, 10])
ax.set_ylim(0, 0.8)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

# --- KDEs
g1 = sns.kdeplot(beta_S_clean, fill=False, color=sns.xkcd_rgb['denim blue'], linestyle='-',
                 label=rf'$\beta_S$ = {med_s:.2f} (10–100 yrs, {n_s} records)', ax=ax)
ax.axvline(x=med_s, ymin=0, ymax=0.1, linewidth=1, color=sns.xkcd_rgb['denim blue'], linestyle='-')

g2 = sns.kdeplot(beta_M_clean, fill=False, color=sns.xkcd_rgb['medium green'], linestyle='-',
                 label=rf'$\beta_M$ = {med_m:.2f} (100–1000 yrs, {n_m} records)', ax=ax)
ax.axvline(x=med_m, ymin=0, ymax=0.1, linewidth=1, color=sns.xkcd_rgb['medium green'], linestyle='-')

g3 = sns.kdeplot(beta_L_clean, fill=False, color=sns.xkcd_rgb['pale red'], linestyle='-',
                 label=rf'$\beta_L$ = {med_l:.2f} (1–10 kyr, {n_l} records)', ax=ax)
ax.axvline(x=med_l, ymin=0, ymax=0.1, linewidth=1, color=sns.xkcd_rgb['pale red'], linestyle='-')

# --- Labels & legend
ax.legend(fontsize=13, loc='upper right', bbox_to_anchor=(1, 1))
ax.set_ylabel('KDE', fontsize=14)
ax.set_xlabel(r'$\beta$', fontsize=14)

# --- Grid & spine styling (dashed grid, black axes)
ax.grid(True, which='major', linestyle='--', linewidth=0.8, color='gray', alpha=0.6)
ax.grid(True, which='minor', linestyle='--', linewidth=0.5, color='gray', alpha=0.4)
for spine in ['bottom', 'left']:
    ax.spines[spine].set_color('black')
    ax.spines[spine].set_linewidth(1.2)

plt.tight_layout()
plt.show()

idx_nan_S = np.where(np.isnan(beta_S))[0]
idx_nan_M = np.where(np.isnan(beta_M))[0]
idx_nan_L = np.where(np.isnan(beta_L))[0]

# Collect labels for NaNs in each array
labels_nan_S = [geoms.series_list[i].label for i in idx_nan_S]
labels_nan_M = [geoms.series_list[i].label for i in idx_nan_M]
labels_nan_L = [geoms.series_list[i].label for i in idx_nan_L]

geoms_S=geoms.copy()
for item in labels_nan_S:
    geoms_S.remove(item)

geoms_M=geoms.copy()
for item in labels_nan_M:
    geoms_M.remove(item)

geoms_L=geoms.copy()
for item in labels_nan_L:
    geoms_L.remove(item)

geoms_S_hue = geoms_S.copy()

for i, ser in enumerate(geoms_S_hue.series_list):
    ser.beta = beta_S_clean[i]

fig,ax = geoms_S_hue.map(hue = 'beta', cmap = 'viridis', projection='Robinson', proj_default={"central_longitude":0})
ax['cb'].set_ylim([0, 5])
ax['cb'].set_ylabel(r'$\beta$', fontsize=14)

geoms_M_hue = geoms_M.copy()

for i, ser in enumerate(geoms_M_hue.series_list):
    ser.beta = beta_M_clean[i]

fig,ax = geoms_M_hue.map(hue = 'beta', cmap = 'viridis', projection='Robinson', proj_default={"central_longitude":0})
ax['cb'].set_ylim([-1, 4])
ax['cb'].set_ylabel(r'$\beta$', fontsize=14)

geoms_L_hue = geoms_L.copy()

for i, ser in enumerate(geoms_L_hue.series_list):
    ser.beta = beta_L_clean[i]

fig,ax = geoms_L_hue.map(hue = 'beta', cmap = 'viridis', projection='Robinson', proj_default={"central_longitude":0})
ax['cb'].set_ylim([-6.5, 6.5])
ax['cb'].set_ylabel(r'$\beta$', fontsize=14)