%% Matlab pseudocode for Retention Clock Matrix (RCM) Analysis
%
% For the accompanying manuscript, please see:
% Defne, Z, NK Ganju, A Aretxabaleta (2016), Estimating time-dependent connectivity in marine systems.
% Geophysical Research Letters
%
% October 2015
% zdefne@usgs.gov
% Zafer Defne
%
--------------------------------------------------------------------------
% START of copying permission statement
%--------------------------------------------------------------------------
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
%--------------------------------------------------------------------------
% END of copying permission statement
%--------------------------------------------------------------------------
clear
%% PART 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Release times shifted to t=0
% read LTRANS output from 24 individual runs
nrelease=24;
% for each release (release label icd)
for icd =1:nrelease
% calculate shift. 2 represents 1 hr shift with 0.5 hr output
shft=2*icd;
shft_last=2*nrelease;
% load lon lat values in LTRANS output
LON{icd+1}=lon( 1+shft : end+1+shft - shft_last , : );
LAT{icd+1}=lat( 1+shft : end+1+shft - shft_last , : );
end
%% PART 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create a stack of connectivity matrices
% read the initial release coordinates
x0=LON(1,:)
y0=LAT(1,:)
% Set number of days to use from LTRANS output
dayn=28;
% number of LTRANS output time steps per day
outPerDay=48;
% total number of time steps
ti_tot=dayn*outPerDay
% LTRANS output time step (1 for 30 mins)
ti_step=1;
% load ESA boundaries
% polyListESA is a list of ESA polygons with x, y vertices for each polygon
% numPolys = total number of Polys
polyListESA{numESAS}(x, y)
% mark the particles in the source polygons initially
for i=1:numPolys
m{1,i}= inpolygon(x0, y0, polyListESA{i}(x), polyListESA{i}(y));
% number of particles in source
S(i)=sum(m{1,i});
end
% mark the particles in the destination polygons at each time step
k=0;
% for each time step
for ti=1 : ti_step : ti_tot
x2=LON(ti,:);
y2=LAT(ti,:);
% mark the particles in the destination polygons
for i=1:numPolys
m{2,i}= inpolygon(x2, y2, polyListESA{i}(x), polyListESA{i}(y));
end
for i=1:numPolys %Source on y axis
for j=1:numPolys %Destination on x axis
% number of particles in destiantion
N(i,j)=sum(m{1,i}.*m{2,j});
N(N==0)=nan;
% connectivity at each time step calculated with eqn: (N_di / N_so)
NbyS(i,j)=N(i,j)/S(i); %
end
end
k=k+1;
% create a stack of connectivity matrices. Pn(time, source, destination)
Pn(k,:,:)=NbyS(:,:);
end
%% PART 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot Retention Clock Matrix (RCM)
% this part of the code is presented mostly as is for ease of reproducibility
% maximum time scale in days
T=28
% time granularity (size of slices in each RCM)
dt=48
%calculate number of timesteps
hrsPerDay=24;
numtstep= T*hrsPerDay / dt;
% calculate averaging window length
wlen=dt*(outPerDay/24)
% calculate average P in each slice
Pavg=[];
for k=1:numtstep
ind=(k-1)*wlen +1 : k*wlen-1;
Pavg(k, :, :)=nanmean(Pn(ind,:,:),1);
end
% set minimum treshold for P (optional)
Pmin=0.001;
Pavg(Pavg0 % print Retention Clock if P is larger than a treshold value
for i=1:numtstep
if(P(i)~=0)
patch([sin(t(i:i+1)) 0], [cos(t(i:i+1)) 0], P(i));
else
patch([sin(t(i:i+1)) 0], [cos(t(i:i+1)) 0], [1 1 1]);
end
end
if di==si
text(0,0, sprintf('%d', si), 'background', 'w', ...
'HorizontalAlignment', 'center', 'fontName', 'Arial')
end
end
set(gca, 'xColor', [1 1 1 ])
set(gca, 'yColor', [1 1 1 ])
set(gca, 'xtick', [])
set(gca, 'ytick', [])
set(gca, 'xtickLabel', [])
set(gca, 'ytickLabel', [])
caxis([0 cmax])
set(gca, 'box', 'on')
if si==di
set(gca,'color', [.5 .5 .5])
else
set(gca,'color', [.8 .8 .8])
end
end
end
% colormap(flipud(copper))
set(gcf, 'paperPositionMode','auto')
set(gcf, 'renderer', 'painters')
% set colormap
colormap(w2c2k(32))
hc = colorbar
set(hc,'position', [ 0.94 0.33 0.015 0.42])
figname=fullfile('connectivity',dirname,sprintf('RCM_%dd_%dhr', T, dt));
% save as Matlab figure
saveas(gcf, figname)
% print eps
print(gcf, '-r600', '-depsc', figname)