@author: Aniket Mitra
PostDoc ISTeP-CYU and AMGC-VUB
---------------------------
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.nonparametric.smoothers_lowess import lowess

# Load the Excel file
file_path = 'filepath'  
sheet_name = "PK02"  
df = pd.read_excel(file_path, sheet_name=sheet_name)

# Ensure column names are stripped of spaces or hidden characters
df.columns = [col.strip() for col in df.columns]

# Define X-axis (distance)
x = df["Distance"]  

# Select trace element columns
selected_elements = [
    "Na/Ca (mmol/mol)", "Mg/Ca (mmol/mol)",
    "Cu/Ca (mmol/mol)", "Zn/Ca (mmol/mol)", "Sr/Ca (mmol/mol)", "Ba/Ca (mmol/mol)", 
    "La/Ca (mmol/mol)", "Ce/Ca (mmol/mol)", "Pr/Ca (mmol/mol)", "Nd/Ca (mmol/mol)", 
    "Sm/Ca (mmol/mol)", "Pb/Ca (mmol/mol)"
]  

# Filter elements that exist in the dataset and reverse the order
elements = [col for col in selected_elements if col in df.columns][::-1]

# Define colors for elements
colors = sns.color_palette("husl", len(elements))

# Create stacked plots (adding 1 for isotopes panel)
fig, axes = plt.subplots(len(elements) + 1, 1, sharex=True, figsize=(8, 18), gridspec_kw={'hspace': 0})

# Plot all selected elements with LOWESS smoothing in reverse order
for i, (element, color) in enumerate(zip(elements, colors)):
    y = df[element].copy()
    

    element_label = element

    # LOWESS smoothing
    smoothed = lowess(y, x, frac=0.1, it=3, return_sorted=False)
    
    ax = axes[i]
    ax.plot(x, y, color=color, alpha=0.2)  # Raw data (transparent)
    ax.plot(x, smoothed, color=color, linewidth=2)  # Smoothed curve
    
    # Scale Y-axis to LOWESS curve range
    ax.set_ylim(np.min(smoothed), np.max(smoothed))  
    
    # Remove gridlines
    ax.yaxis.grid(False)
    ax.xaxis.grid(False)

    # Formatting
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.spines["left"].set_visible(False)
    
    # Alternating Y-axis labels
    if i % 2 == 0:
        ax.yaxis.set_label_position("left")
        ax.yaxis.tick_left()
    else:
        ax.yaxis.set_label_position("right")
        ax.yaxis.tick_right()
    
    ax.set_ylabel(f"{element_label}", rotation=90, labelpad=15, fontsize=12, color=color, va='center')
    
    # Add vertical lines at tick marks and minor ticks
    ax.yaxis.set_minor_locator(plt.AutoLocator())
    ax.tick_params(axis='y', which='both', length=4, direction='in', right=True)
    ax.yaxis.set_major_formatter(plt.ScalarFormatter())
    ax.yaxis.set_minor_formatter(plt.ScalarFormatter())
    ax.tick_params(axis='y', which='minor', length=2, direction='in', right=True)

# Add δ¹³C and δ¹⁸O in the final panel with separate Y-axes and box plots
ax_iso = axes[-1]
ax_iso2 = ax_iso.twinx()

# Use 'Distance stable' column for isotopes
x_stable = df["Distance stable"]

# Group data by Distance stable and compute statistics
d13C_stats = df.groupby("Distance stable")["d13C (‰)"].describe()
d18O_stats = df.groupby("Distance stable")["d18O (‰)"].describe()

# Boxplot for δ¹³C
d13C_positions = d13C_stats.index
d13C_values = [df[df["Distance stable"] == pos]["d13C (‰)"].values for pos in d13C_positions]
ax_iso.boxplot(d13C_values, positions=d13C_positions, widths=0.1, patch_artist=True, boxprops=dict(facecolor="blue", alpha=0.3))
ax_iso.plot(d13C_stats.index, d13C_stats["mean"], color='blue', label="δ¹³C (‰)")

# Boxplot for δ¹⁸O
d18O_positions = d18O_stats.index
d18O_values = [df[df["Distance stable"] == pos]["d18O (‰)"].values for pos in d18O_positions]
ax_iso2.boxplot(d18O_values, positions=d18O_positions, widths=0.1, patch_artist=True, boxprops=dict(facecolor="red", alpha=0.3))
ax_iso2.plot(d18O_stats.index, d18O_stats["mean"], color='red', label="δ¹⁸O (‰)")

# Reverse Y-axis for both
ax_iso.set_ylim(ax_iso.get_ylim()[::-1])  # Reverse Y-axis for δ¹³C
ax_iso2.set_ylim(ax_iso2.get_ylim()[::-1])  # Reverse Y-axis for δ¹⁸O

ax_iso.set_ylabel("δ¹³C (‰)", fontsize=12, color='blue')
ax_iso2.set_ylabel("δ¹⁸O (‰)", fontsize=12, color='red')


# Add minor tick marks for Y-axis in isotope panel
ax_iso.yaxis.set_minor_locator(plt.AutoLocator())
ax_iso.tick_params(axis='y', which='minor', length=2, direction='in', left=True)
ax_iso2.yaxis.set_minor_locator(plt.AutoLocator())
ax_iso2.tick_params(axis='y', which='minor', length=2, direction='in', right=True)

# X-axis label only at the bottom
axes[-1].set_xlabel("Distance (mm)")

plt.show()
