% Script to produce figure 6 from Hernandez et al 2020
% author: Laia Comas-Bru
% ORCID: https://orcid.org/0000-0002-7882-4996 
% Matlab version 2019b
% Date last edited: June 2020
% script available at 10.5281/zenodo.3898382 
% necessary files: 
% - sea-level pressure dataset from 20CRv2c (Compo et al., 2011):
% "prmsl.csv"
% - brewermap function for Matlab: https://github.com/DrosteEffect/BrewerMap 
% - "Sites_map.m" script (if uncommented)

clear
close all
%%
X=importdata('prmsl.csv',',');Coord=X(:,1:2); A=X(:,3:1682);

%Rearrange dataset(delete months between Apr-Nov incl.)
mSLP=ones(2556,417);
for k=1:139
    clear winter
    winter=ones(2556,3);
    n=12*k;
      
for i=1:3 % 3 because it is DJF
winter(:,i)=A(:,n+i-1);
end
mSLP(:,3*k-2:3*k)=winter(:,:);
end

mSLP=mSLP/100;% Convert wSLP from Pa to hPa

% Length of the record Des1871-Feb2009: winters 1872-2009
mSLP(:,415:417)=[];% 2556*417(414=3*138yrs)

wSLP=ones(2556,138);
clear i
for i=1:138 % number of years
    clear a;
    a= mSLP(:,1+(3*(i-1)):3+(3*(i-1)));
    wSLP(:,i)=mean(a,2);
end

%% load NAO and EA indices

% load NAO winter composite (DJF; Comas-Bru and Hernandez, 2018)
opts = delimitedTextImportOptions("NumVariables", 2);
opts.DataLines = [2, Inf];
opts.Delimiter = ",";
opts.VariableNames = ["year", "NAO"];
opts.VariableTypes = ["double", "double"];
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
NAO = readtable("[ADD_YOUR_PATH_HERE]\NAO_winter_composite.csv", opts);
clear opts

% load EA winter composite (DJF; Comas-Bru and Hernandez, 2018)
opts = delimitedTextImportOptions("NumVariables", 2);
opts.DataLines = [2, Inf];
opts.Delimiter = ",";
opts.VariableNames = ["year", "EA"];
opts.VariableTypes = ["double", "double"];
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
EA = readtable("[ADD_YOUR_PATH_HERE]\EA_winter_composite.csv", opts);
clear opts

EA(1:50,:)=[];NAO(1:50,:)=[];% years not in CRU removed: Period 1902-2016
EA=table2array(EA); NAO=table2array(NAO);
clear a;a(:,1)=NAO(:,2);a(:,2)=EA(:,2);

% NAO/EA COMBINATION OF SIGNS
clear m;[m]=find(a(:,1)>0);    Pos=m; % NAO>0    
clear m;[m]=find(a(:,1)>0 & a(:,2)>0);    PosPos=m; % NAO>0 and EA>0
clear m;[m]=find(a(:,1)>0 & a(:,2)<0);    PosNeg=m;% NAO>0 and EA<0
clear m;[m]=find(a(:,1)<0);    Neg=m; % NAO<0
clear m;[m]=find(a(:,1)<0 & a(:,2)>0);    NegPos=m; % NAO<0 and EA>0
clear m;[m]=find(a(:,1)<0 & a(:,2)<0);    NegNeg=m; % NAO<0 and EA<0  

Psame=length(PosPos);Nsame=length(NegNeg);
same=PosPos; same(Psame+1:Nsame+Psame,1)=NegNeg(:,1);

Pop=length(PosNeg);Nop=length(NegPos);
op=PosNeg; opp(Pop+1:Nop+Pop)=NegPos(:,1);

clear Psame Pop m Pos* Neg* Nil* Nsame Nop

%% Figures
load mymap_grey
load coastlines; 

% Fig. 3,1 SAME
clear i
samewSLP=ones(2556,length(same));
for i=1:length(same)
clear k; k=same(i);
samewSLP(:,i)=wSLP(:,k);
end

clear a c z; a=mean(samewSLP,1); c=size(samewSLP,2); z=ones(size(samewSLP));

for i=1:c
   z(:,i)=samewSLP(:,i)-a(1,i);
end

CorrMatrix=corr(z'); c=min(CorrMatrix,[],2); % (2556,1)
c=abs(c); %clear n; [n]=find(c>0.85);  % to find the locations of CAs
a1=(1:1:71)';Coord(:,4)=repmat(a1,36,1); a2=(10:2:80)';a2=flipud(a2);

clear i
for i=1:length(a2)
[m]=find(Coord(:,2)==a2(i));
Coord(m,3)=i; % Latitude
end
clear cMatrix
cMatrix=ones(36,71);
clear i
for i=1:2556
cMatrix(Coord(i,3),Coord(i,4))=c(i,1);
end
cMatrix=flipud(cMatrix);%clear m;[m]=find(cMatrix<0.6);cMatrix(m)=NaN*ones(size(m));

% FIND MAXs

% clear i
% for i=1:length(n)
% LabLat(i,1)=Coord(n(i),2);
% LabLon(i,1)=Coord(n(i),1);
% end
% 
% [nj]=find(LabLon(:,1)>300);
% clear i
% for i=1:length(nj)
%     k=nj(i);
%     LabLon(k,1)=LabLon(k,1)-360;
% end

% [Ice]=find(LabLat>50);LabelLatIce=mean(LabLat(Ice,1));LabelLonIce=mean(LabLon(Ice,1));
% [Azo]=find(LabLat<50);LabelLatAzo=mean(LabLat(Azo,1));LabelLonAzo=mean(LabLon(Azo,1));

latlim=[10 80];lonlim=[-100 40]; 

figure; subplot(3,1,1); 
axesm('lambertstd','MapLatLimit',latlim,'MapLonLimit',lonlim,...
    'Frame','on','Grid','on','MeridianLabel','off','ParallelLabel','off')
axis off; 
setm(gca,'MLabelLocation',60,'FLineWidth',1);
surfm(latlim,lonlim,cMatrix);
contourm(latlim,lonlim,cMatrix,0.6:0.1:0.9,'k');
% textm(LabelLatIce, LabelLonIce, '+','Color','r','Fontsize',20);
% textm(LabelLatAzo, LabelLonAzo, '+','Color','r','Fontsize',20);
t = title('Same');set(t, 'FontSize', 10,'FontWeight','bold');
plotm(coastlat,coastlon,'k')
set(gcf,'renderer','Painters')
caxis([0.5,1]); colorbar; colormap(mymap_grey(2:end-2,:))
% Sites_map


% Fig. 3,2 ALL DATASET
clear a c z; a=mean(wSLP,1); c=size(wSLP,2); z=ones(size(wSLP));

for i=1:c
   z(:,i)=wSLP(:,i)-a(1,i);
end

CorrMatrix=corr(z'); c=min(CorrMatrix,[],2); % (2556,1)
c=abs(c); clear n; [n]=find(c>0.85);  % to find the locations of CAs

% Create matrix with Coord
a1=(1:1:71)';Coord(:,4)=repmat(a1,36,1); a2=(10:2:80)';a2=flipud(a2);

clear i
for i=1:length(a2)
[m]=find(Coord(:,2)==a2(i));
Coord(m,3)=i; % Latitude
end
clear cMatrix
cMatrix=ones(36,71);
clear i
for i=1:2556
cMatrix(Coord(i,3),Coord(i,4))=c(i,1);
end
cMatrix=flipud(cMatrix);%clear m;[m]=find(cMatrix<0.6);cMatrix(m)=NaN*ones(size(m));

% FIND MAXs

% clear i
% for i=1:length(n)
% LabLat(i,1)=Coord(n(i),2);
% LabLon(i,1)=Coord(n(i),1);
% end
% 
% [nj]=find(LabLon(:,1)>300);
% clear i
% for i=1:length(nj)
%     k=nj(i);
%     LabLon(k,1)=LabLon(k,1)-360;
% end

% [Ice]=find(LabLat>50);LabelLatIce=mean(LabLat(Ice,1));LabelLonIce=mean(LabLon(Ice,1));
% [Azo]=find(LabLat<50);LabelLatAzo=mean(LabLat(Azo,1));LabelLonAzo=mean(LabLon(Azo,1));

gcf; subplot(3,1,2); 
axesm('lambertstd','MapLatLimit',latlim,'MapLonLimit',lonlim,...
    'Frame','on','Grid','on','MeridianLabel','off','ParallelLabel','off')
axis off; 
setm(gca,'MLabelLocation',60,'FLineWidth',1);
surfm(latlim,lonlim,cMatrix);
contourm(latlim,lonlim,cMatrix,0.6:0.1:0.9,'k');
% textm(LabelLatIce, LabelLonIce, '+','Color','r','Fontsize',20);
% textm(LabelLatAzo, LabelLonAzo, '+','Color','r','Fontsize',20);
t = title('All');set(t, 'FontSize', 10,'FontWeight','bold');
plotm(coastlat,coastlon,'k')
set(gcf,'renderer','Painters')
caxis([0.5,1]); colorbar; colormap(mymap_grey(2:end-2,:))
% Sites_map

%% Fig. 3,3 OP
clear i
opSLP=ones(2556,length(op));
for i=1:length(op)
clear k; k=op(i);
opSLP(:,i)=wSLP(:,k);
end

clear a c z; a=mean(opSLP,1); c=size(opSLP,2); z=ones(size(opSLP));

for i=1:c
   z(:,i)=opSLP(:,i)-a(1,i);
end

CorrMatrix=corr(z'); c=min(CorrMatrix,[],2); % (2556,1)
c=abs(c); clear n; [n]=find(c>0.85);  % to find the locations of CAs

% Create matrix with Coord
a1=(1:1:71)';Coord(:,4)=repmat(a1,36,1); a2=(10:2:80)';a2=flipud(a2);

clear i
for i=1:length(a2)
[m]=find(Coord(:,2)==a2(i));
Coord(m,3)=i; % Latitude
end
clear cMatrix
cMatrix=ones(36,71);
clear i
for i=1:2556
cMatrix(Coord(i,3),Coord(i,4))=c(i,1);
end
cMatrix=flipud(cMatrix);%clear m;[m]=find(cMatrix<0.6);cMatrix(m)=NaN*ones(size(m));

% FIND MAXs

% clear i
% for i=1:length(n)
% LabLat(i,1)=Coord(n(i),2);
% LabLon(i,1)=Coord(n(i),1);
% end
% 
% [nj]=find(LabLon(:,1)>300);
% clear i
% for i=1:length(nj)
%     k=nj(i);
%     LabLon(k,1)=LabLon(k,1)-360;
% end
% 
% [Ice]=find(LabLat>50);LabelLatIce=mean(LabLat(Ice,1));LabelLonIce=mean(LabLon(Ice,1));
% [Azo]=find(LabLat<50);LabelLatAzo=mean(LabLat(Azo,1));LabelLonAzo=mean(LabLon(Azo,1));

gcf; subplot(3,1,3); 
axesm('lambertstd','MapLatLimit',latlim,'MapLonLimit',lonlim,...
    'Frame','on','Grid','on','MeridianLabel','off','ParallelLabel','off')
axis off; 
setm(gca,'MLabelLocation',60,'FLineWidth',1);
surfm(latlim,lonlim,cMatrix);
[c,h] = contourm(latlim,lonlim,cMatrix,0.6:0.1:0.9,'k');
% textm(LabelLatIce, LabelLonIce, '+','Color','r','Fontsize',20);
% textm(LabelLatAzo, LabelLonAzo, '+','Color','r','Fontsize',20);
t = title('op');set(t, 'FontSize', 10,'FontWeight','bold');
plotm(coastlat,coastlon,'k')
set(gcf,'renderer','Painters')
caxis([0.5,1]); colorbar; colormap(mymap_grey(2:end-2,:))
% Sites_map

%% save figure
print('-bestfit',sprintf('Figure_6.pdf'),'-dpdf')

clear
close all
