Drawing Pollen Diagrams with MATLAB – MATLAB and Python Recipes for Earth Sciences

The 1st edition of MDRES includes a MATLAB script called multiplot to create diagrams in a design typically used by palynologists and micropaleontologists. Such diagrams to display several pollen or microfossil records use patches together with horizontal lines or bars on axes with the same scales. I modified the script using the new string array and the new text property rotation, both introduced in R2016b, to make the graphics look even more like a pollen diagram. The new script will be included in the 2nd edition of MDRES to be submitted to the publisher before the end of October 2017.

The history of the development of my script goes back a few years, when a doctoral student asked about the software Tilia, developed by Eric Grimm, Director of Sciences Emeritus of the Illinois State Museum. However, the student did not want to use the full catalogue of functions of the excellent software, but only to display their diatom counts, and therefore I hesitated to purchase the software. We first need to create some synthetic data. The first column contains the age vector, while all other data are uniformly-distributed pseudo random numbers. The total of all values within each column is 100. We can change the number of data points m.

clear, clc, close all
rng(10)
m = 40;
data(:,1) = sortrows(m*rand(m,1));
data(:,2:15) = rand(m,14); 
data(:,2) = 5*data(:,2);
data(:,3) = 5*data(:,3);
data(:,4) = 10*data(:,4);
data(:,5) = 5*data(:,5);
data(:,6) = 5*data(:,6);
data(:,7) = 3*data(:,7);
data(:,8) = 5*data(:,8);
data(:,9) = 5*data(:,9);
data(:,10) = 15*data(:,10);
data(:,11) = 12*data(:,11);
data(:,12) = 5*data(:,12);
data(:,13) = 3*data(:,13);
data(:,14) = 2*data(:,14);
data(:,15) = 2*data(:,15);
data(1,1) = 0;

We also create a list of fossil names by typing

names = ["Lorem ipsum"
 "Dolor sit"
 "Amet consectetuer"
 "Adipiscing elit"
 "Aenean commodo"
 "Ligula eget"
 "Dolor aenean"
 "Massa cum"
 "Sociis natoque"
 "Penatibus et"
 "Magnis dis"
 "Parturient montes"
 "Nascetur ridiculus"
 "Mus donec"
 "Quam felis"
];

Alternatively, we can import any other data set with an arbitrary number of data points and variables, provided the total comes to 100%. We add an age of zero at both ends of the time series by typing

[m,n] = size(data);
data(2:m+1,:) = data;
data(1,:) = zeros(1,n);
data(1,1) = data(2,1);
data(m+2,:) = zeros(1,n);
data(m+2,1) = data(m+1,1);

The multiplot function first creates a figure window with the same size as the screen. It then plots all records at the same scale.

root1 = groot;
scrsz = root1.ScreenSize;
scrsz = 0.6*scrsz;
sf = 0.75;
colorcode = [0.8 0.8 0.4];
figure1 = figure(...
    'Position',[1 scrsz(4) scrsz(3) scrsz(4)],...
    'Color',[1 1 1]);
for i = 1:n-1
if i == 1
xcoord = max(data(:,2))/100; 
axes1 = axes(...
    'Parent',figure1,...
    'YMinorTick','on',...
    'XMinorTick','on',...
    'TickDir','out',...
    'Position',[0.1 0.1 sf*xcoord 0.6],...
    'FontSize',12,...
    'XLim',[0 max(data(:,2))],...
    'YLim',[min(data(:,1)) max(data(:,1))],...
    'YDir','reverse',...
    'Box','off');
patch1 = patch(data(:,2),data(:,1),colorcode);
hold on
%line1 = line(data(:,2),data(:,1),...
%    'Color',[0 0 0]);
barh1 = barh(data(2:end-1,1),data(2:end-1,2),...
    'BarWidth',0,...
    'EdgeColor',[0 0 0],...
    'FaceColor',colorcode);
title1 = title(names(i),...
    'Rotation',45,...
    'FontWeight','normal',...
    'Position',[0 -1],...
    'HorizontalAlignment','left',...
    'VerticalAlignment','top',...
    'FontSize',12);
xlabel1 = ylabel('Age (cal yrs BP)');
ylabel1 = xlabel('%');
xcoord = xcoord + 0.01;
else 
axes2 = axes(...
    'Parent',figure1,...
    'YTick',zeros(1,0),...
    'FontSize',12,...
    'YMinorTick','on',...
    'XMinorTick','on',...
    'TickDir','out',...
    'Position',[sf*xcoord+0.1 0.1 ...
    sf*max(data(:,i+1))/100 0.6],...
    'XLim',[0 max(data(:,i+1))],...
    'YLim',[min(data(:,1)) max(data(:,1))],...
    'YDir','reverse',...
    'Box','off');
 patch2 = patch(data(:,i+1),data(:,1),colorcode);
 hold on
 %line2 = line(data(2:end-1,i+1),data(2:end-1,1),...
 %    'Color',[0 0 0]);
 barh2 = barh(data(2:end-1,1),data(2:end-1,i+1),...
    'BarWidth',0,...
    'EdgeColor',[0 0 0],...
    'FaceColor',colorcode);
 xlabel2 = xlabel('%');
 title(names(i),...
    'Rotation',45,...
    'FontWeight','normal',...
    'Position',[0 -1],...
    'HorizontalAlignment','left',...
    'VerticalAlignment','top',...
    'FontSize',12);
 xcoord = xcoord + max(data(:,i+1))/100 + 0.01;
end
end

We can change the multiplier for scrsz between 0 and 1 to match the figure size with the screen size. We can also change colorcode to modify the color of the patches. The variable sf scales the graph with respect to the figure window. Instead of a patch plot, the graphics function line for a line plot is also available. To use either of these we need to uncomment the patch or the line command by removing the % signs at the beginning of the corresponding lines.