% 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: 
% - temperature and precipitation data from CRU (Harris et al., 2014):
% "cru_ts4.01.1901.2016.pre.dat.nc" and "cru_ts4.01.1901.2016.tmp.dat.nc"
% - brewermap function for Matlab: https://github.com/DrosteEffect/BrewerMap 
% - "Sites_map.m" script (if uncommented)

clear
close all

%% load CRU data

clear path;
path=('[ADD_YOUR_PATH_HERE]\cru_ts4.01.1901.2016.pre.dat.nc');
pre=ncread(path,'pre');%mm/month
time=ncread(path,'time'); time=double(time);%months*1901-2016
lon=ncread(path,'lon'); lon=double(lon); lat=ncread(path,'lat');lat=double(lat);

clear path; path=('[ADD_YOUR_PATH_HERE]\cru_ts4.01.1901.2016.tmp.dat.nc');
tmp=ncread(path,'tmp');%degC

clear path
Tbase = datenum(1900, 1, 1);
duration=sscanf(datestr(time(end,1) + Tbase),'16-Dec-%d')-sscanf(datestr(time(1,1) + Tbase),'16-Jan-%d')+1; %years

time_mm(:,1)=repmat((1:1:12)',duration,1);
time_mm(:,2)=sort(repmat((sscanf(datestr(time(1,1) + Tbase),'16-Jan-%d'):sscanf(datestr(time(end,1) + Tbase),'16-Dec-%d'))',12,1),'ascend');

clear Tbase time

%empty cells to NaN
clear m; [m]=find(pre==9.969209968386869e+36);
if isempty(m)==0
    pre(m)=NaN;
end

clear m; [m]=find(tmp==9.969209968386869e+36);
if isempty(m)==0
    tmp(m)=NaN;
end

clear m;

% Region: latlim=[10 80];lonlim=[-100 40];
latlim=[30 75];lonlim=[-100 50];

clear m; [m]=find(lon>min(lonlim) & lon<max(lonlim));clear n; [n]=find(lat>min(latlim) & lat<max(latlim));
tmp=tmp(m,n,:); pre=pre(m,n,:);lat=lat(n,1); lon=lon(m,1);

%1. delete months Mar to Nov (3-11) => DJF
clear m; [m]=find(time_mm(:,1)~=12 & time_mm(:,1)~=1 & time_mm(:,1)~=2);
pre(:,:,m)=[]; tmp(:,:,m)=[]; time_mm(m,:)=[];

%2. delete months at the edge
pre(:,:,1:2)=[]; tmp(:,:,1:2)=[]; time_mm(1:2,:)=[];
pre(:,:,end)=[]; tmp(:,:,end)=[]; time_mm(end,:)=[];

% calculate DJFM averages
clear pre_res; pre_res=reshape(pre, size(pre,1)*size(pre,2),size(pre,3));
clear tmp_res; tmp_res=reshape(tmp, size(tmp,1)*size(tmp,2),size(tmp,3));

clear x; x=ones(size(pre_res,1),(size(time_mm,1)/3))*NaN;
clear y; y=ones(size(tmp_res,1),(size(time_mm,1)/3))*NaN;

for st=1:size(pre_res,1) %each gridcell
    for k=1:(size(time_mm,1)/3) %each year
        x(st,k)=nansum(pre_res(st,3*k-2:3*k))/3; % winter precip in mm/month
        y(st,k)=nanmean(tmp_res(st,3*k-2:3*k)); % winter temp
    end
end

TWP=reshape(x,size(pre,1),size(pre,2),size(x,2)); % Total Winter Precip
MWT=reshape(y,size(tmp,1),size(tmp,2),size(y,2)); % Mean Winter Temperature

TWP(TWP==0)=NaN;

clear x y st k

time=(min(time_mm(:,2))+1:max(time_mm(:,2)))';%new time

%% ----
% 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);
opp=PosNeg; opp(Pop+1:Nop+Pop)=NegPos(:,1);

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

%% correlations
% latlim1=[30 70];lonlim1=[-100 50];
 load coastlines;
 cmapw(1:1024,:)=brewermap(1024,'RdYlBu');  
            
%NAO/EA OPPOSITE
clear A B1* B2*
A=NAO(opp,:);
B1a=MWT(:,:,opp); B1=reshape(B1a, size(B1a,1)*size(B1a,2),size(B1a,3)); clear B1a
B2a=TWP(:,:,opp);B2=reshape(B2a, size(B2a,1)*size(B2a,2),size(B2a,3));clear B2a

length_opp=size(opp,1)

RHO=ones(size(B1,1),1)*NaN; PVALUES=ones(size(B1,1),1)*NaN;
for i=1:size(B1,1)
clear X;
X(:,1)=A(:,2);
X(:,2)=B1(i,:);
[rho,pval]=corr(X,'type','Spearman'); 
RHO(i,1)=rho(2,1);
PVALUES(i,1)=pval(2,1);
end

clear m; [m]=find(PVALUES>0.05); RHO(m)=NaN;
NAO_TMP_OPP=reshape(RHO,size(MWT,1), size(MWT,2));

            figure; subplot(2,3,2);
            set(gcf,'renderer','Painters')
            axesm('lambertstd','MapLatLimit',latlim,'MapLonLimit',lonlim);
            axis off; framem on;setm(gca,'FLineWidth',1);gridm on;mlabel off;plabel off;
            surfm(latlim,lonlim,flipud(rot90(NAO_TMP_OPP)));
            plotm(coastlat,coastlon,'k')
            title('NAO TMP OPP'); colorbar; caxis([-1,1]);colormap(flipud(cmapw))
%             Sites_map

RHO=ones(size(B2,1),1)*NaN; PVALUES=ones(size(B2,1),1)*NaN;
for i=1:size(B2,1)
clear X;
X(:,1)=A(:,2);
X(:,2)=B2(i,:);
[rho,pval]=corr(X,'type','Spearman'); 
RHO(i,1)=rho(2,1);
PVALUES(i,1)=pval(2,1);
end
clear m; [m]=find(PVALUES>0.05); RHO(m)=NaN;
NAO_PRE_OPP=reshape(RHO,size(TWP,1), size(TWP,2));

            gcf; subplot(2,3,5);
            set(gcf,'renderer','Painters')
            axesm('lambertstd','MapLatLimit',latlim,'MapLonLimit',lonlim);
            axis off; framem on;setm(gca,'FLineWidth',1);gridm on;mlabel off;plabel off;
            surfm(latlim,lonlim,flipud(rot90(NAO_PRE_OPP)));
            plotm(coastlat,coastlon,'k')
            title('NAO PRE OPP'); colorbar; caxis([-1,1]);colormap(flipud(cmapw))
%             Sites_map

%NAO/EA SAME
clear A B1* B2*
A=NAO(same,:);
B1a=MWT(:,:,same); B1=reshape(B1a, size(B1a,1)*size(B1a,2),size(B1a,3)); clear B1a
B2a=TWP(:,:,same);B2=reshape(B2a, size(B2a,1)*size(B2a,2),size(B2a,3));clear B2a

length_same=size(same,1)

RHO=ones(size(B1,1),1)*NaN; PVALUES=ones(size(B1,1),1)*NaN;
for i=1:size(B1,1)
clear X;
X(:,1)=A(:,2);
X(:,2)=B1(i,:);
[rho,pval]=corr(X,'type','Spearman'); 
RHO(i,1)=rho(2,1);
PVALUES(i,1)=pval(2,1);
end
clear m; [m]=find(PVALUES>0.05); RHO(m)=NaN;
NAO_TMP_SAME=reshape(RHO,size(MWT,1), size(MWT,2));

            gcf; subplot(2,3,3); 
            set(gcf,'renderer','Painters')
            axesm('lambertstd','MapLatLimit',latlim,'MapLonLimit',lonlim);
            axis off; framem on;setm(gca,'FLineWidth',1);gridm on;mlabel off;plabel off;
            surfm(latlim,lonlim,flipud(rot90(NAO_TMP_SAME)));
            plotm(coastlat,coastlon,'k')
            title('NAO TMP SAME'); colorbar; caxis([-1,1]);colormap(flipud(cmapw))
%             Sites_map
            
RHO=ones(size(B2,1),1)*NaN; PVALUES=ones(size(B2,1),1)*NaN;
for i=1:size(B2,1)
clear X;
X(:,1)=A(:,2);
X(:,2)=B2(i,:);
[rho,pval]=corr(X,'type','Spearman'); 
RHO(i,1)=rho(2,1);
PVALUES(i,1)=pval(2,1);
end
clear m; [m]=find(PVALUES>0.05); RHO(m)=NaN;
NAO_PRE_SAME=reshape(RHO,size(TWP,1), size(TWP,2));

            gcf; subplot(2,3,6); 
            set(gcf,'renderer','Painters')
            axesm('lambertstd','MapLatLimit',latlim,'MapLonLimit',lonlim);
            axis off; framem on;setm(gca,'FLineWidth',1);gridm on;mlabel off;plabel off;
            surfm(latlim,lonlim,flipud(rot90(NAO_PRE_SAME)));
            plotm(coastlat,coastlon,'k')
            title('NAO PRE SAME'); colorbar; caxis([-1,1]);colormap(flipud(cmapw))
%             Sites_map
            
%NAO/EA ALL
clear A B1* B2*
A=NAO;
B1a=MWT;B1=reshape(B1a, size(B1a,1)*size(B1a,2),size(B1a,3)); clear B1a
B2a=TWP;B2=reshape(B2a, size(B2a,1)*size(B2a,2),size(B2a,3)); clear B2a

RHO=ones(size(B1,1),1)*NaN; PVALUES=ones(size(B1,1),1)*NaN;
for i=1:size(B1,1)
clear X;
X(:,1)=A(:,2);
X(:,2)=B1(i,:);
[rho,pval]=corr(X,'type','Spearman'); 
RHO(i,1)=rho(2,1);
PVALUES(i,1)=pval(2,1);
end
clear m; [m]=find(PVALUES>0.05); RHO(m)=NaN;
NAO_TMP_ALL=reshape(RHO,size(MWT,1), size(MWT,2));

            gcf; subplot(2,3,1); 
            set(gcf,'renderer','Painters')
            axesm('lambertstd','MapLatLimit',latlim,'MapLonLimit',lonlim);
            axis off; framem on;setm(gca,'FLineWidth',1);gridm on;mlabel off;plabel off;
            surfm(latlim,lonlim,flipud(rot90(NAO_TMP_ALL)));
            plotm(coastlat,coastlon,'k')
            title('NAO TMP ALL'); colorbar; caxis([-1,1]);colormap(flipud(cmapw))
%             Sites_map
            
RHO=ones(size(B2,1),1)*NaN; PVALUES=ones(size(B2,1),1)*NaN;
for i=1:size(B2,1)
clear X;
X(:,1)=A(:,2);
X(:,2)=B2(i,:);
[rho,pval]=corr(X,'type','Spearman'); 
RHO(i,1)=rho(2,1);
PVALUES(i,1)=pval(2,1);
end
clear m; [m]=find(PVALUES>0.05); RHO(m)=NaN;
NAO_PRE_ALL=reshape(RHO,size(TWP,1), size(TWP,2));

            gcf; subplot(2,3,4); 
            set(gcf,'renderer','Painters')
            axesm('lambertstd','MapLatLimit',latlim,'MapLonLimit',lonlim);
            axis off; framem on;setm(gca,'FLineWidth',1);gridm on;mlabel off;plabel off;
            surfm(latlim,lonlim,flipud(rot90(NAO_PRE_ALL)));
            plotm(coastlat,coastlon,'k')
            title('NAO PRE ALL'); colorbar; caxis([-1,1]);colormap(flipud(cmapw))
%             Sites_map

%% save figure
h=gcf; set(h,'PaperOrientation','landscape');
print('-bestfit',sprintf('Figure_7.pdf'),'-dpdf')

clear
close all
            