Written by: Kaustubh Thirumalai, Brown University (Github)
Written on: November 21, 2018
Make sure the archived KL-126 dataset is in your working folder prior to running the code.
import matplotlib.pyplot as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
df = pd.read_table('KL-126-Pangea.txt').values
age_d18O = df[:,2]
d18O = df[:,3]
A line plot of the raw record appears to replicate Figure 2b in the Kudrass et al. 2001 Geology paper.
sns.set(style='darkgrid',palette='muted',font_scale=1.5)
[fig,ax] = plt.subplots(1,1,figsize=(11,5))
p = ax.plot(age_d18O,d18O,color='xkcd:mid blue')
ax.invert_yaxis()
ax.set_title('KL-126 Record')
ax.set_ylabel('$G. ruber$ δ$\mathregular{^{18}}$O (‰)')
ax.set_xlabel('Age (ka)')
plt.tight_layout()
plt.show()
#---- Save figure ----#
fig.savefig('Fig1.pdf')
Now, the problem becomes more apparent if we plot this record using a step function. For this we can use the drawstyle='steps-mid'
option in matplotlib. Finally, the 'interpolation' problem becomes apparent when we plot markers at each data point in the archived record. You can see swaths of datapoints that are clearly interpolated from one measurement to another.
[fig,(ax1,ax2)] = plt.subplots(2,1,figsize=(11,10))
p1 = ax1.plot(age_d18O,d18O,'-',color='xkcd:mid blue',drawstyle='steps-mid')
ax1.invert_yaxis()
ax1.set_title('KL-126 Record')
ax1.set_ylabel('$G. ruber$ δ$\mathregular{^{18}}$O (‰)')
ax1.set_xlabel('Age (ka)')
p2 = ax2.plot(age_d18O,d18O,'o-',color='xkcd:mid blue',drawstyle='steps-mid',markerfacecolor='xkcd:light orange')
ax2.invert_yaxis()
ax2.set_ylabel('$G. ruber$ δ$\mathregular{^{18}}$O (‰)')
ax2.set_xlabel('Age (ka)')
plt.tight_layout()
plt.show()
#---- Save figure ----#
fig.savefig('Fig2.pdf')
Let's look at the Last Glacial Maximum to present (20,000 yrs ago to modern)
[fig,(ax1,ax2)] = plt.subplots(2,1,figsize=(11,10))
p1 = ax1.plot(age_d18O,d18O,'-',color='xkcd:mid blue',drawstyle='steps-mid')
ax1.invert_yaxis()
ax1.set_title('KL-126 Record')
ax1.set_ylabel('$G. ruber$ δ$\mathregular{^{18}}$O (‰)')
ax1.set_xlabel('Age (ka)')
ax1.set_xlim([-1,20])
p2 = ax2.plot(age_d18O,d18O,'o-',color='xkcd:mid blue',drawstyle='steps-mid',markerfacecolor='xkcd:light orange')
ax2.invert_yaxis()
ax2.set_ylabel('$G. ruber$ δ$\mathregular{^{18}}$O (‰)')
ax2.set_xlabel('Age (ka)')
ax2.set_xlim([-1,20])
plt.tight_layout()
plt.show()
#---- Save figure ----#
fig.savefig('Fig3.pdf')
Something seems off!
We can plot the difference of each timestep (first difference function) to see what the resolution of the archived data is and how it changes over time.
first_diff = np.diff(age_d18O)*1000
first_age = age_d18O[1:]
[fig,(ax1,ax2)] = plt.subplots(2,1,figsize=(11,10))
# Plot of archived record
p1 = ax1.plot(age_d18O,d18O,'-',color='xkcd:mid blue',drawstyle='steps-mid')
ax1.invert_yaxis()
ax1.set_title('KL-126 Record')
ax1.set_ylabel('$G. ruber$ δ$\mathregular{^{18}}$O (‰)')
ax1.set_xlabel('Age (ka)')
# Plot of resolution over time
p2 = ax2.plot(first_age,first_diff,'-',color='xkcd:purple',drawstyle='steps-mid')
p3 = ax2.plot(first_age,50*np.ones(len(first_age)),'--',color='xkcd:red',drawstyle='steps-mid')
ax2.annotate('50 years',[0,20],color='xkcd:red')
ax2.set_title('Time resolution of archived data')
ax2.set_ylabel('Time step (years)')
ax2.set_xlabel('Age (ka)')
plt.tight_layout()
plt.show()
#---- Save figure ----#
fig.savefig('Fig4.pdf')
We see that it changes over time but it is also plateaued at times! We can now observe how the original timeseries was interpolated -- to the point that the record sometimes has sub-decadal resolution. The age model precludes such incredible sedimentation rates and certainly, something is awry.
We can approximate the "original" measurements by finding local minima and maxima in the archived (interpolated) record. Of course, these are not going to mimic the actual measurments one-to-one, but unfortunately, since the actual measurements are not archived, this will have to do.
from scipy.signal import find_peaks
peaks, _ = find_peaks(d18O)
peaksy, _ = find_peaks(-1*d18O)
[fig,ax1] = plt.subplots(1,1,figsize=(15,5))
# Plot of archived record
p1 = ax1.plot(age_d18O,d18O,'-',color='xkcd:mid blue')
ax1.invert_yaxis()
ax1.set_title('KL-126 Record')
ax1.set_ylabel('$G. ruber$ δ$\mathregular{^{18}}$O (‰)')
ax1.set_xlabel('Age (ka)')
ax1.plot(age_d18O[peaks],d18O[peaks],'o',color='xkcd:bright red',markersize=10,markerfacecolor='none')
ax1.plot(age_d18O[peaksy],d18O[peaksy],'o',color='xkcd:green',markersize=10,markerfacecolor='none')
ax1.set_xlim([-1,105])
plt.tight_layout()
plt.show()
#---- Save figure ----#
fig.savefig('Fig5.pdf')
Now we will put in all the green and red points from above into one array and retain their indices in time to get the approximated measurements.
# Join the arrays
p = np.transpose((age_d18O[peaks],d18O[peaks]))
n = np.transpose((age_d18O[peaksy],d18O[peaksy]))
K = np.concatenate((p,n),axis=0)
# Add the first point
K = np.transpose([np.append(K[:,0],age_d18O[0]),np.append(K[:,1],d18O[0])])
# Sort by time after joining arrays
L = K[K[:,0].argsort()]
mod_age = L[:,0]
mod_d18O = L[:,1]
[fig,ax1] = plt.subplots(1,1,figsize=(11,5))
# Plot of archived record
p1 = ax1.plot(age_d18O,d18O,'-',color='xkcd:mid blue',linewidth=2.5,label='Archived Data')
ax1.invert_yaxis()
ax1.set_title('KL-126 Record')
ax1.set_ylabel('$G. ruber$ δ$\mathregular{^{18}}$O (‰)')
ax1.set_xlabel('Age (ka)')
ax1.plot(mod_age,mod_d18O,'o-',color='xkcd:orange',label='Reconstructed Data')
ax1.legend()
plt.tight_layout()
plt.show()
#---- Save figure ----#
fig.savefig('Fig6.pdf')
It's not perfect - but it's better than the interpolated dataset in my opinion!
Now let's focus on the alkenone dataset: first, let's convert the data from being on depth to age. We can do this by interpolating the age model.
depth = df[:,0]
age_depth = df[:,1]
alk_depth = df[:,4]
alk_SST = df[:,5]
alk_age = np.interp(alk_depth,depth,age_depth)
[fig,(ax1,ax2)] = plt.subplots(2,1,figsize=(8,8))
p1 = ax1.plot(alk_depth,alk_SST,'-s',color='xkcd:purple blue',markersize=4)
ax1.set_ylabel('Alkenone SST (°C)')
ax1.set_xlabel('Depth (m)')
p2 = ax2.plot(alk_age,alk_SST,'-s',color='xkcd:purple',markersize=4)
ax2.set_ylabel('Alkenone SST (°C)')
ax2.set_xlabel('Age (ka)')
plt.tight_layout()
plt.show()
#---- Save figure ----#
fig.savefig('Fig7.pdf')
[fig,(ax1,ax2)] = plt.subplots(2,1,figsize=(12,8))
Xax = [-1,101]
p1 = ax1.plot(mod_age,mod_d18O,'o-',color='xkcd:orange',label="δ$\mathregular{^{18}}$O",markersize=5,linewidth=1)
ax1.set_ylabel('$G. ruber$ δ$\mathregular{^{18}}$O')
ax1.set_xlim(Xax)
ax1.set_xlabel('Age (ka)')
ax1.xaxis.set_label_position('top')
ax1.invert_yaxis()
ax1.tick_params(
axis='x', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
bottom=False, # ticks along the bottom edge are off
top=False, # ticks along the top edge are off
labeltop= True,
labelbottom=False) # labels along the bottom edge are off
ax1.grid(axis='y',linestyle='--')
ax1.grid(axis='x',linestyle='--')
# Sea-Level
pos1 = ax1.get_position() # get the original position
pos2 = [pos1.x0,pos1.y0-(pos1.height),pos1.width,pos1.height+0.005]
ax2.set_position(pos2) # set a new position
p2 = ax2.plot(alk_age,alk_SST,'s-',color='xkcd:purple',label="Sea Level",markersize=4)
ax2.set_xlim(Xax)
ax2.set_ylabel('SST (°C)')
ax2.set_xlabel('Age (ka)')
ax2.grid(axis='y',linestyle='--')
ax2.grid(axis='x',linestyle='--')
plt.show()
#---- Save figure ----#
fig.savefig('Fig8.pdf')