Displaying Sediment Records with both Age and Depth Axis with MATLAB

Sediment data is typically display in diagrams with two axes, a core depth axis and an age axis. Here is a MATLAB script to convert depths into ages and display a sediment parameter in a diagram with both age and depth axis.

First, we create a synthetic data set of a sine wave with a period of 200.

clear, clc, close all
t = 0 : 1000; 
rec = 0.5*sin(2*pi*t/200) + 0.5;

Then we load the age depth conversion table contained in the file tiepoints.txt.

tiepoints = load('tiepoints.txt');
tiepoints = 1.7*tiepoints;

We can display the relationship between age and depth in an age-depth plot by typing

figure('Position',[100 800 500 500],...
  'Color',[1 1 1])
axes('Position',[0.2 0.2 0.7 0.7],...
  'XLim',[0 1000],...
  'YLim',[0 500],...
  'YDir','Reverse',...
  'Box','On',...
  'FontSize',18)
line(tiepoints(:,2),tiepoints(:,1),...
  'LineWidth',1)
xlabel('Age (kyr BP)')
ylabel('Depth in Core (m)')

Next, we interpolate this relationship to an evenly-spaced age scale.

age = 0 : 20 : 1000; age = age';
depthint = interp1(tiepoints(:,2),tiepoints(:,1),...
   age,'pchip');
depthintlabels = num2str(depthint,'%3.0f\n');

and display in a first diagram with two axes.

figure1 = figure('Position',[100 500 1800 600],...
  'Color',[1 1 1]) 
ax(1) = axes('Position',[0.1 0.4 0.8 0.4],...
  'Color','None',...
  'XLim',[0 1000],...
  'YLim',[-0.1 1.1],...
  'XTick',0:50:1000,...
  'YLim',[-0.1 1.1],...
  'FontSize',18)
line1 = line(t,rec,...
  'LineWidth',1);
xlabel(ax(1),'Age (kyr BP)')
ylabel(ax(1),'Proxy Value')
ax(2) = axes('Position',[0.1 0.25 0.8 0.4],...
  'Color','None',...
  'XLim',[0 1000],...
  'XTickMode','Manual',...
  'XTick',0:50:1000,...
  'XTickLabels',depthintlabels,...
  'YLim',[0 1],...
  'YTick',[],...
  'YColor','None',...
  'FontSize',18)
xlabel(ax(2),'Depth in Core (m)')

Then we interpolate the  ages to an evenly-spaced depth scale

depth = 0 : 20 : 500;
ageint = interp1(depthint,age,depth,'pchip');
depthlabels = num2str(depth,'%3.0f\n');

and again display in a first diagram with two axes.

figure1 = figure('Position',[100 100 1800 600],...
  'Color',[1 1 1]);
ax(1) = axes('Position',[0.1 0.4 0.8 0.4],...
  'Color','None',...
  'XLim',[0 1000],...
  'XTick',0:50:1000,...
  'YLim',[-0.1 1.1],...
  'FontSize',18);
line1 = line(t,rec,...
  'LineWidth',1);
xlabel(ax(1),'Age (kyr BP)')
ylabel(ax(1),'Proxy Value')
ax(2) = axes('Position',[0.1 0.25 0.8 0.4],...
  'Color','None',...
  'XLim',[0 1000],...
  'XTickMode','Manual',...
  'XTick',ageint,...
  'XTickLabels',depthlabels,...
  'YLim',[0 1],...
  'YTick',[],...
  'YColor','None',...
  'FontSize',18);
xlabel(ax(2),'Depth in Core (m)')