#!/usr/bin/env python
# coding: utf-8

# # The curious case of KL-126
# 
# #### *Program to resample the archived KL-126 δ¹⁸O dataset in order to simulate the "actual measurements"*
# 
# ----
# Written by: Kaustubh Thirumalai, Brown University ([Github](https://github.com/holy-kau))  
# Written on: November 21, 2018  
# 
# [Make sure the archived KL-126 dataset is in your working folder prior to running the code.](https://kthi.squarespace.com/s/KL-126-PANGEA.txt)

# In[1]:


import matplotlib.pyplot as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns


# ### Import KL-126 δ¹⁸O dataset and load into Pandas dataframe

# In[2]:


df = pd.read_table('KL-126-Pangea.txt').values
age_d18O = df[:,2]
d18O = df[:,3]


# ### Line plot of raw record: So what's the problem?
# 
# A line plot of the raw record appears to replicate Figure 2b in the Kudrass et al. 2001 Geology paper.

# In[3]:


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')


# ### A closer look
# 
# 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.

# In[4]:


[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')


# ### An even closer look in time!
# Let's look at the Last Glacial Maximum to present (20,000 yrs ago to modern)

# In[5]:


[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!

# ### What's the resolution?
# 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. 

# In[6]:


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.

# ### Approximating the original measurements
# 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.

# In[7]:


from scipy.signal import find_peaks


# In[8]:


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')


# ### Array manipulation
# 
# 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.

# In[9]:


# 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!

# ## Alkenones
# 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.

# In[10]:


depth = df[:,0]
age_depth = df[:,1]
alk_depth = df[:,4]
alk_SST = df[:,5]

alk_age = np.interp(alk_depth,depth,age_depth)


# In[11]:


[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')


# ## Final Plot

# In[12]:


[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')