
%*************************************************************
%*   ENSO Impacts Child Undernutrition in the Global Tropics
%*	Jesse K. Anttila-Hughes, Amir S. Jina, Gordon C. McCord
%*	Corresponding author: Amir Jina, amirjina@uchicago.edu
%*   This .m file produces the pixel-level teleconnection values 
%*  that define the sample in the paper, and specifically figures
%*  1c and 1d
%*************************************************************/	

clear
%load in udel data

cd 'data directory'

data = importdata('tas_Amon_UDel_historical_halfdeg_1950-2014_DATA.txt');

% Data should be formatted in long-format with columns:
% longitude
% latitude
% month (with January 1950 as 1)
% temperature

%NOTE: the code is identical for both temperature and precipitation, but is performed after importing
% the pr_Amon_UDel_historical_halfdeg_1950-2014_DATA.txt file instead of the tas_* file above

data = data.data;
year = floor(data(:,3)/12+1950);
month = mod(data(:,3),12);
month(month==0)=12;

%run corrs for each pixel
%total vector length
total_length = length([min(data(:,2)):0.5:max(data(:,2))])*length([min(data(:,1)):0.5:max(data(:,1))]);

%make monthly correlation matrices
ind_M = find(month==1);
corr_data = data(ind_M,:);
month_array=NaN(total_length,12,2+2014-1949);
test_lat = reshape(corr_data(:,2),total_length,2014-1949);
test_lon = reshape(corr_data(:,1),total_length,2014-1949);
month_array(:,:,1)=test_lon(:,1:12);
month_array(:,:,2)=test_lat(:,1:12);

for iM = 1:12
    ind_M = find(month==iM);
    corr_data = data(ind_M,:);
    month_array(:,iM,3:end) = reshape(corr_data(:,4),total_length,1,2014-1949);
end


%load in and format nino3.4 data
% data obtained from ERSSTv5 (centered base periods): http://www.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/detrend.nino34.ascii.txt
%data are minimally formatted to make them readable by matlab

nino34 = importdata('nino34_formatted.csv');
nino34 = nino34(81:145,:);
nino34_corr = NaN(total_length,12,2014-1949);

for iM = 1:12
n_temp=permute(nino34(:,iM+1),[3 2 1]);
nino34_corr(:,iM,:)=repmat(n_temp,total_length,1,1);
end 

%format two-month lag
cd 'C:\Users\Amir\Dropbox\enso_anthro_working_dir\climate_data\nino'
nino34_o = importdata('nino34_formatted.csv');
nino34_L2(:,1)=nino34_o(2:end,1);
for iR=2:length(nino34_o)
    nino34_L2(iR-1,2:13)=[nino34_o(iR-1,12:13) nino34_o(iR,2:11)];
end
nino34_L2 = nino34_L2(80:144,:);
nino34_corr_L2 = NaN(total_length,12,2014-1949);
for iM = 1:12
n_temp=permute(nino34_L2(:,iM+1),[3 2 1]);
nino34_corr_L2(:,iM,:)=repmat(n_temp,total_length,1,1);
end 

%correct for that longitude misalignment (-180 to 180 instead of 0 to 360)
month_array_temp = month_array(:,1,1);
month_array_m360 = month_array(:,1,1)-360;
month_array_recentered = NaN(length(month_array(:,1,1)),1);

for i=1:length(month_array(:,1,1))
    if month_array_temp(i,1)>180
        month_array_recentered(i,1)=month_array_m360(i,1);
    else
        month_array_recentered(i,1)=month_array_temp(i,1);
    end
end

% calculate correlation
%loop
pix_corr=NaN(total_length,14);
pix_corr(:,1)=month_array_recentered;
pix_corr(:,2)=month_array(:,1,2);

pix_pval=NaN(total_length,14);
pix_pval(:,1)=month_array_recentered;
pix_pval(:,2)=month_array(:,1,2);

pix_corr_L2=NaN(total_length,14);
pix_corr_L2(:,1)=month_array_recentered;
pix_corr_L2(:,2)=month_array(:,1,2);

pix_pval_L2=NaN(total_length,14);
pix_pval_L2(:,1)=month_array_recentered;
pix_pval_L2(:,2)=month_array(:,1,2);


for iM=1:12
    iM
   for iX=1:total_length
      [R,P]=corrcoef(month_array(iX,iM,3:end),nino34_corr(iX,iM,:)); 
      pix_corr(iX,iM+2)=R(2,1);
      pix_pval(iX,iM+2)=P(2,1);
      
      [R,P]=corrcoef(month_array(iX,iM,3:end),nino34_corr_L2(iX,iM,:)); 
      pix_corr_L2(iX,iM+2)=R(2,1);
      pix_pval_L2(iX,iM+2)=P(2,1);
   end
end

pix_ssig=pix_pval(:,3:end);
pix_ssig(pix_ssig<0.1)=1;
pix_ssig(pix_ssig~=1)=0;
pix_ssig_tot=sum(pix_ssig,2);
pix_ssig=[pix_pval(:,1:2) pix_ssig];
pix_ssig_tot=[pix_pval(:,1:2) pix_ssig_tot];

pix_ssig_L2=pix_pval_L2(:,3:end);
pix_ssig_L2(pix_ssig_L2<0.1)=1;
pix_ssig_L2(pix_ssig_L2~=1)=0;
pix_ssig_tot_L2=sum(pix_ssig_L2,2);
pix_ssig_L2=[pix_pval_L2(:,1:2) pix_ssig_L2];
pix_ssig_tot_L2=[pix_pval_L2(:,1:2) pix_ssig_tot_L2];

%calculate sign of correlation for L2
pix_ssig_corr_L2 = [pix_pval(:,1:2) pix_ssig_L2(:,3:14).*pix_corr_L2(:,3:14)];
pix_ssig_corr_L2(pix_ssig_corr_L2==0)=NaN;
pix_ssig_corr_L2(:,3:14)=sign(pix_ssig_corr_L2(:,3:14));
pix_ssig_posneg_L2 = [pix_pval(:,1:2) zeros(length(pix_pval),2)];
for i=1:length(pix_pval)
    pix_ssig_posneg_L2(i,3) = sum(pix_ssig_corr_L2(i,:)==1);
    pix_ssig_posneg_L2(i,4) = sum(pix_ssig_corr_L2(i,:)==-1);
end
temp1 = pix_ssig_posneg_L2(:,3);
temp2 = pix_ssig_posneg_L2(:,4); 
temp1(temp1<3)=0;
temp1(temp1>=3)=1;
temp2(temp2<3)=0;
temp2(temp2>=3)=1;
pix_posneg_dummy_L2 = [pix_pval(:,1:2) temp1 temp2];

%make csv for mapping and definition of country teleconnections
cd 'output directory'

cHeader = {'lon' 'lat' 'pos_corr' 'neg_corr'}; %dummy header
commaHeader = [cHeader;repmat({','},1,numel(cHeader))]; %insert commaas
commaHeader = commaHeader(:)';
textHeader = cell2mat(commaHeader); %cHeader in text with commas
%write header to file
fid = fopen('nino34_udel_posnegdummy_L2.csv','w'); 
fprintf(fid,'%s\n',textHeader)
fclose(fid)
%write data to end of file
dlmwrite('nino34_udel_posnegdummy_L2.csv',pix_posneg_dummy_L2,'-append');

cHeader = {'lon' 'lat' 'jan_pvl' 'feb_pvl' 'mar_pvl' 'apr_pvl' 'may_pvl' 'jun_pvl' 'jul_pvl' 'aug_pvl' 'sep_pvl' 'oct_pvl' 'nov_pvl' 'dec_pvl'}; %dummy header
commaHeader = [cHeader;repmat({','},1,numel(cHeader))]; %insert commaas
commaHeader = commaHeader(:)';
textHeader = cell2mat(commaHeader); %cHeader in text with commas
%write header to file
fid = fopen('nino34_udel_ssig_mth.csv','w'); 
fprintf(fid,'%s\n',textHeader)
fclose(fid)
%write data to end of file
dlmwrite('nino34_udel_ssig_mth.csv',pix_ssig,'-append');

cHeader = {'lon' 'lat' 'jan_pvl' 'feb_pvl' 'mar_pvl' 'apr_pvl' 'may_pvl' 'jun_pvl' 'jul_pvl' 'aug_pvl' 'sep_pvl' 'oct_pvl' 'nov_pvl' 'dec_pvl'}; %dummy header
commaHeader = [cHeader;repmat({','},1,numel(cHeader))]; %insert commaas
commaHeader = commaHeader(:)';
textHeader = cell2mat(commaHeader); %cHeader in text with commas
%write header to file
fid = fopen('nino34_udel_ssig_pvl.csv','w'); 
fprintf(fid,'%s\n',textHeader)
fclose(fid)
%write data to end of file
dlmwrite('nino34_udel_ssig_pvl.csv',pix_pval,'-append');

cHeader = {'lon' 'lat' 'jan_pvl' 'feb_pvl' 'mar_pvl' 'apr_pvl' 'may_pvl' 'jun_pvl' 'jul_pvl' 'aug_pvl' 'sep_pvl' 'oct_pvl' 'nov_pvl' 'dec_pvl'}; %dummy header
commaHeader = [cHeader;repmat({','},1,numel(cHeader))]; %insert commaas
commaHeader = commaHeader(:)';
textHeader = cell2mat(commaHeader); %cHeader in text with commas
%write header to file
fid = fopen('nino34_udel_ssig_cor.csv','w'); 
fprintf(fid,'%s\n',textHeader)
fclose(fid)
%write data to end of file
dlmwrite('nino34_udel_ssig_cor.csv',pix_corr,'-append');

cHeader = {'lon' 'lat' 'jan_pvl' 'feb_pvl' 'mar_pvl' 'apr_pvl' 'may_pvl' 'jun_pvl' 'jul_pvl' 'aug_pvl' 'sep_pvl' 'oct_pvl' 'nov_pvl' 'dec_pvl'}; %dummy header
commaHeader = [cHeader;repmat({','},1,numel(cHeader))]; %insert commaas
commaHeader = commaHeader(:)';
textHeader = cell2mat(commaHeader); %cHeader in text with commas
%write header to file
fid = fopen('nino34_udel_ssig_mth_L2.csv','w'); 
fprintf(fid,'%s\n',textHeader)
fclose(fid)
%write data to end of file
dlmwrite('nino34_udel_ssig_mth_L2.csv',pix_ssig_L2,'-append');

cHeader = {'lon' 'lat' 'jan_pvl' 'feb_pvl' 'mar_pvl' 'apr_pvl' 'may_pvl' 'jun_pvl' 'jul_pvl' 'aug_pvl' 'sep_pvl' 'oct_pvl' 'nov_pvl' 'dec_pvl'}; %dummy header
commaHeader = [cHeader;repmat({','},1,numel(cHeader))]; %insert commaas
commaHeader = commaHeader(:)';
textHeader = cell2mat(commaHeader); %cHeader in text with commas
%write header to file
fid = fopen('nino34_udel_ssig_pvl_L2.csv','w'); 
fprintf(fid,'%s\n',textHeader)
fclose(fid)
%write data to end of file
dlmwrite('nino34_udel_ssig_pvl_L2.csv',pix_pval_L2,'-append');

cHeader = {'lon' 'lat' 'jan_pvl' 'feb_pvl' 'mar_pvl' 'apr_pvl' 'may_pvl' 'jun_pvl' 'jul_pvl' 'aug_pvl' 'sep_pvl' 'oct_pvl' 'nov_pvl' 'dec_pvl'}; %dummy header
commaHeader = [cHeader;repmat({','},1,numel(cHeader))]; %insert commaas
commaHeader = commaHeader(:)';
textHeader = cell2mat(commaHeader); %cHeader in text with commas
%write header to file
fid = fopen('nino34_udel_ssig_cor_L2.csv','w'); 
fprintf(fid,'%s\n',textHeader)
fclose(fid)
%write data to end of file
dlmwrite('nino34_udel_ssig_cor_L2.csv',pix_corr_L2,'-append');

cHeader = {'lon' 'lat' 'sigcor_tot'}; %dummy header
commaHeader = [cHeader;repmat({','},1,numel(cHeader))]; %insert commaas
commaHeader = commaHeader(:)';
textHeader = cell2mat(commaHeader); %cHeader in text with commas
%write header to file
fid = fopen('nino34_udel_ssig_tot.csv','w'); 
fprintf(fid,'%s\n',textHeader)
fclose(fid)
%write data to end of file
dlmwrite('nino34_udel_ssig_tot.csv',pix_ssig_tot,'-append');

cHeader = {'lon' 'lat' 'sigcor_tot'}; %dummy header
commaHeader = [cHeader;repmat({','},1,numel(cHeader))]; %insert commaas
commaHeader = commaHeader(:)';
textHeader = cell2mat(commaHeader); %cHeader in text with commas
%write header to file
fid = fopen('nino34_udel_ssig_tot_L2.csv','w'); 
fprintf(fid,'%s\n',textHeader)
fclose(fid)
%write data to end of file
dlmwrite('nino34_udel_ssig_tot_L2.csv',pix_ssig_tot_L2,'-append');

