% Code to compute the final steps in the Vegetation Sensitivity Index

% FILE DESCRIPTION

% anom+"var"+.csv = anomalies from mean-variance regression for a given variable
% anomtemp.csv
% anomaet.csv
% anomcld.csv
% anomevi.csv

% coef+"var"+.csv = coefficients from regression from a given variable
% coeftemp.csv
% coefaet.csv
% coefcld.csv
% coeft1.csv lag t-1 coefficient

% coef+NonT1_global_+"var"+.csv = coefficients from regression for a given variable -
% T1 not accounted in the weight - Globally weighted - for RGB only
% coefNonT1_global_temp.csv
% coefNonT1_global_aet.csv
% coefNonT1_global_cld.csv

% coef+NonT1+"var"+.csv = coefficients from regression for a given variable -
% T1 not accounted in the weight - Locally weighted - for Individual
% display
% coefNonT1temp.csv
% coefNonT1aet.csv
% coefNonT1cld.csv
% 
% coef+"Perc"+"var"+.csv = percent coefficient from regression for a given variable
% coefPerct1.csv
% coefPerctemp.csv
% coefPercaet.csv
% coefPerccld.csv
%
% Number of months per pixel used in the regression
% sigMonth.csv
%
% cI + var.csv Multiple regression confidence intervals for each variable
% cIcld.csv
% cIaet.csv
% cItemp.csv
% cIt1.csv
%
% the raw coefficients - in order to divide cIs by them to scale them
% cI_coefcld.csv
% cI_coefaet.csv
% cI_coeftemp.csv
% cI_coeft1.csv
%
% Mean of the monthly standard errors (SD/mean) of the number of valid
% MODIS observations
% month_yr_stderror.csv

% ----------------------------------------------------------------------%

% 1. Import variability files normalised by mean vs. variance
EVI=csvread('anomevi.csv');EVI(EVI==-9999)=NaN;EVI(EVI==-8888)=NaN;
Cld=csvread('anomcld.csv');Cld(Cld==-9999)=NaN;Cld(Cld==-8888)=NaN;
Tmp=csvread('anomtemp.csv');Tmp(Tmp==-9999)=NaN;Tmp(Tmp==-8888)=NaN;
AET=csvread('anomaet.csv');AET(AET==-9999)=NaN;AET(AET==-8888)=NaN;

% 2. Compute non-weighted Sensitivity (not weighted according to the coefficients)
% 2.1.Compute individual non-weighted sensitivities
SensCldNW=log10((EVI+1)./(Cld+1));SensCldNW=(SensCldNW-min(min(SensCldNW)))./(max(max(SensCldNW))-min(min(SensCldNW))).*100;
SensTmpNW=log10((EVI+1)./(Tmp+1));SensTmpNW=(SensTmpNW-min(min(SensTmpNW)))./(max(max(SensTmpNW))-min(min(SensTmpNW))).*100;
SensAETNW=log10((EVI+1)./(AET+1));SensAETNW=(SensAETNW-min(min(SensAETNW)))./(max(max(SensAETNW))-min(min(SensAETNW))).*100;
% 2.2 Compute Overall non-weighted Sensitivity
a=reshape(SensCldNW,3000*7200,1);
b=reshape(SensTmpNW,3000*7200,1);
c=reshape(SensAETNW,3000*7200,1);
M=[a b c];
SensTotalNW=nansum(M,2); clear M;
a(a>-1)=1;
b(b>-1)=1;
c(c>-1)=1;
M=[a b c];clear a; clear b;clear c;
M=nansum(M,2);
M(M>0)=1;
M(M==0)=NaN;
SensTotalNW=SensTotalNW.*M;
SensTotalNW=reshape(SensTotalNW,3000,7200);
SensTotalNW=(SensTotalNW-min(min(SensTotalNW)))./(max(max(SensTotalNW))-min(min(SensTotalNW))).*100;

% 3. Import the NT1 coefficient files - globally weighted - for RGB & Weighted Sensitivity
CoeffCldG=csvread('coefNonT1_global_cld.csv');CoeffCldG(CoeffCldG==-9999)=NaN;CoeffCldG(CoeffCldG==-8888)=NaN;CoeffCldG=reshape(CoeffCldG,3000*7200,1);CoeffCldG=CoeffCldG.*M;CoeffCldG=reshape(CoeffCldG,3000,7200);
CoeffAETG=csvread('coefNonT1_global_aet.csv');CoeffAETG(CoeffAETG==-9999)=NaN;CoeffAETG(CoeffAETG==-8888)=NaN;CoeffAETG=reshape(CoeffAETG,3000*7200,1);CoeffAETG=CoeffAETG.*M;CoeffAETG=reshape(CoeffAETG,3000,7200);
CoeffTmpG=csvread('coefNonT1_global_temp.csv');CoeffTmpG(CoeffTmpG==-9999)=NaN;CoeffTmpG(CoeffTmpG==-8888)=NaN;CoeffTmpG=reshape(CoeffTmpG,3000*7200,1);CoeffTmpG=CoeffTmpG.*M;CoeffTmpG=reshape(CoeffTmpG,3000,7200);
% 3.1 Re-scaling of the coefficients to 0-1 in order to use them in the
% weighted Sensitivity
CoeffCldGW=CoeffCldG./100;
CoeffAETGW=CoeffAETG./100;
CoeffTmpGW=CoeffTmpG./100;

% 4. Compute coefficient-weighted Sensitivity (weighted according to regression coefficients)
% 4.1.Compute individual coefficient-weighted sensitivities
SensCldW=SensCldNW.*CoeffCldGW;
SensTmpW=SensTmpNW.*CoeffTmpGW;
SensAETW=SensAETNW.*CoeffAETGW;
% 4.2 Compute overall coefficient-weighted Sensitivity (i.e. VSI)
a=reshape(SensCldW,3000*7200,1);
b=reshape(SensTmpW,3000*7200,1);
c=reshape(SensAETW,3000*7200,1);
W=[a b c];clear a b c
SensTotalW=nansum(W,2);clear W;
SensTotalW=SensTotalW.*M;
SensTotalW=reshape(SensTotalW,3000,7200);
SensTotalW=(SensTotalW-min(min(SensTotalW)))./(max(max(SensTotalW))-min(min(SensTotalW))).*100;
% 4.3 Rescale coefficient-weighted sensitivities to 0 -100 for display
% purposes
SensCldW=(SensCldW-min(min(SensCldW)))./(max(max(SensCldW))-min(min(SensCldW))).*100;
SensTmpW=(SensTmpW-min(min(SensTmpW)))./(max(max(SensTmpW))-min(min(SensTmpW))).*100;
SensAETW=(SensAETW-min(min(SensAETW)))./(max(max(SensAETW))-min(min(SensAETW))).*100;

% 5. Import the NT1 coefficient files - individual variable weighted - for individual display
CoeffCldL=csvread('coefNonT1cld.csv');CoeffCldL(CoeffCldL==-9999)=NaN;CoeffCldL(CoeffCldL==-8888)=NaN;CoeffCldL=reshape(CoeffCldL,3000*7200,1);CoeffCldL=CoeffCldL.*M;CoeffCldL=reshape(CoeffCldL,3000,7200);
CoeffTmpL=csvread('coefNonT1temp.csv');CoeffTmpL(CoeffTmpL==-9999)=NaN;CoeffTmpL(CoeffTmpL==-8888)=NaN;CoeffTmpL=reshape(CoeffTmpL,3000*7200,1);CoeffTmpL=CoeffTmpL.*M;CoeffTmpL=reshape(CoeffTmpL,3000,7200);
CoeffAETL=csvread('coefNonT1aet.csv');CoeffAETL(CoeffAETL==-9999)=NaN;CoeffAETL(CoeffAETL==-8888)=NaN;CoeffAETL=reshape(CoeffAETL,3000*7200,1);CoeffAETL=CoeffAETL.*M;CoeffAETL=reshape(CoeffAETL,3000,7200);

% 6. Import the Percent coefficient files - they give an estimate of the
% relative weight of each variable in the regression - including T1
CoeffPercTmp=csvread('coefPerctemp.csv');CoeffPercTmp(CoeffPercTmp==-9999)=NaN;CoeffPercTmp(CoeffPercTmp==-8888)=NaN;CoeffPercTmp=reshape(CoeffPercTmp,3000*7200,1);CoeffPercTmp=CoeffPercTmp.*M;CoeffPercTmp=reshape(CoeffPercTmp,3000,7200);
CoeffPercAET=csvread('coefPercaet.csv');CoeffPercAET(CoeffPercAET==-9999)=NaN;CoeffPercAET(CoeffPercAET==-8888)=NaN;CoeffPercAET=reshape(CoeffPercAET,3000*7200,1);CoeffPercAET=CoeffPercAET.*M;CoeffPercAET=reshape(CoeffPercAET,3000,7200);
CoeffPercCld=csvread('coefPerccld.csv');CoeffPercCld(CoeffPercCld==-9999)=NaN;CoeffPercCld(CoeffPercCld==-8888)=NaN;CoeffPercCld=reshape(CoeffPercCld,3000*7200,1);CoeffPercCld=CoeffPercCld.*M;CoeffPercCld=reshape(CoeffPercCld,3000,7200);
CoeffPercT1=csvread('coefPerct1.csv');CoeffPercT1(CoeffPercT1==-9999)=NaN;CoeffPercT1(CoeffPercT1==-8888)=NaN;CoeffPercT1=reshape(CoeffPercT1,3000*7200,1);CoeffPercT1=CoeffPercT1.*M;CoeffPercT1=reshape(CoeffPercT1,3000,7200);

% 7. Import the file with Number of Significant Months
sigMonth=csvread('sigMonth.csv');sigMonth(sigMonth==-9999)=NaN;sigMonth=reshape(sigMonth,3000*7200,1);sigMonth=sigMonth.*M;sigMonth=reshape(sigMonth,3000,7200);

% 8. Import the Confidence Interval files
cIcld=csvread('cIcld.csv');cIcld(cIcld==-9999)=NaN;cIcld(cIcld==-8888)=NaN;
cIaet=csvread('cIaet.csv');cIaet(cIaet==-9999)=NaN;cIaet(cIaet==-8888)=NaN;
cItmp=csvread('cItemp.csv');cItmp(cItmp==-9999)=NaN;cItmp(cItmp==-8888)=NaN;
cIt1=csvread('cIt1.csv');cIt1(cIt1==-9999)=NaN;cIt1(cIt1==-8888)=NaN;
% 8.1 Import the raw coefficient values - to scale the cI
cI_coefcld=csvread('cI_coefcld.csv');cI_coefcld(cI_coefcld==-9999)=NaN;cI_coefcld(cI_coefcld==-8888)=NaN;
cI_coefaet=csvread('cI_coefaet.csv');cI_coefaet(cI_coefaet==-9999)=NaN;cI_coefaet(cI_coefaet==-8888)=NaN;
cI_coeftmp=csvread('cI_coeftemp.csv');cI_coeftmp(cI_coeftmp==-9999)=NaN;cI_coeftmp(cI_coeftmp==-8888)=NaN;
cI_coeft1=csvread('cI_coeft1.csv');cI_coeft1(cI_coeft1==-9999)=NaN;cI_coeft1(cI_coeft1==-8888)=NaN;
% 8.2 Confidence Intervals tmp
% 8.2.1. Remove extreme values from the cI files (the top 0.01%)
cItmpv=reshape(cItmp,3000*7200,1);
cI_coeftmpv=reshape(cI_coeftmp,3000*7200,1);
p=prctile(cItmpv,99.99);
z=prctile(cI_coeftmpv,99.99);
cItmp(cItmp>=p)=NaN;
cI_coeftmp(cI_coeftmp>=z)=NaN;
% 8.2.2. Scale according to the raw coefficient value
cItmp_S=cItmp./cI_coeftmp;
cItmp_Sv=reshape(cItmp_S,3000*7200,1);
% 8.2.3. Remove the top 1% of the scaled values
x=prctile(cItmp_Sv,99);
cItmp_S(cItmp_S>=x)=NaN;
% 8.3. Confidence Intervals cld
% 8.3.1 Remove extreme values from the cI files (the top 0.01%)
cIcldv=reshape(cIcld,3000*7200,1);
cI_coefcldv=reshape(cI_coefcld,3000*7200,1);
p=prctile(cIcldv,99.99);
z=prctile(cI_coefcldv,99.99);
cIcld(cIcld>=p)=NaN;
cI_coefcld(cI_coefcld>=z)=NaN;
% 8.3.2. Scale according to the raw coefficient value
cIcld_S=cIcld./cI_coefcld;
cIcld_Sv=reshape(cIcld_S,3000*7200,1);
% 8.3.3. Remove the top 1% of the scaled values
x=prctile(cIcld_Sv,99);
cIcld_S(cIcld_S>=x)=NaN;
% 8.4. Confidence Intervals aet
% 8.4.1. Remove extreme values from the cI files (the top 0.01%)
cIaetv=reshape(cIaet,3000*7200,1);
cI_coefaetv=reshape(cI_coefaet,3000*7200,1);
p=prctile(cIaetv,99.99);
z=prctile(cI_coefaetv,99.99);
cIaet(cIaet>=p)=NaN;
cI_coefaet(cI_coefaet>=z)=NaN;
% 8.4.2. Scale according to the raw coefficient value
cIaet_S=cIaet./cI_coefaet;
cIaet_Sv=reshape(cIaet_S,3000*7200,1);
% 8.4.3. Remove the top 1% of the scaled values
x=prctile(cIaet_Sv,99);
cIaet_S(cIaet_S>=x)=NaN;
% 8.5. Confidence Intervals t1
% 8.5.1. Remove extreme values from the cI files (the top 0.01%)
cIt1v=reshape(cIt1,3000*7200,1);
cI_coeft1v=reshape(cI_coeft1,3000*7200,1);
p=prctile(cIt1v,99.99);
z=prctile(cI_coeft1v,99.99);
cIt1(cIt1>=p)=NaN;
cI_coeft1(cI_coeft1>=z)=NaN;
% 8.5.2. Scale according to the raw coefficient value
cIt1_S=cIt1./cI_coeft1;
cIt1_Sv=reshape(cIt1_S,3000*7200,1);
% 8.5.3. Remove the top 1% of the scaled values
x=prctile(cIt1_Sv,99);
cIt1_S(cIt1_S>=x)=NaN;
clear p; clear z; clear x


% 9. Import the SE of the observations
month_yr_stderror=csvread('month_yr_stderror.csv');month_yr_stderror(month_yr_stderror==-9999)=NaN;month_yr_stderror(month_yr_stderror==-8888)=NaN;
% 9.1. Remove extreme values (the top 0.1%)
month_yr_stderror=reshape(month_yr_stderror,3000*7200,1);
month_yr_stderror=month_yr_stderror.*M;
p=prctile(month_yr_stderror,99.9);
month_yr_stderror(month_yr_stderror>p)=NaN;
month_yr_stderror=reshape(month_yr_stderror,3000,7200);
clear p

save Sensitivity

% 10. Export files to ASCII format

% NOTE: if converted to raster in ArcMap, add the following header to the the text files: 
% ncols 7200
% nrows 3000
% xllcorner -180
% yllcorner -60
% cellsize 0.05

SensCldNW(isnan(SensCldNW))=-9999;dlmwrite('SensCldNW.txt',SensCldNW,'\t');
SensTmpNW(isnan(SensTmpNW))=-9999;dlmwrite('SensTmpNW.txt',SensTmpNW,'\t');
SensAETNW(isnan(SensAETNW))=-9999;dlmwrite('SensAETNW.txt',SensAETNW,'\t');
SensTotalNW(isnan(SensTotalNW))=-9999;dlmwrite('SensTotalNW.txt',SensTotalNW,'\t');

SensCldW(isnan(SensCldW))=-9999;dlmwrite('SensCldW.txt',SensCldW,'\t');
SensTmpW(isnan(SensTmpW))=-9999;dlmwrite('SensTmpW.txt',SensTmpW,'\t');
SensAETW(isnan(SensAETW))=-9999;dlmwrite('SensAETW.txt',SensAETW,'\t');
SensTotalW(isnan(SensTotalW))=-9999;dlmwrite('SensTotalW.txt',SensTotalW,'\t');

CoeffCldL(isnan(CoeffCldL))=-9999;dlmwrite('CoeffNT1CldL.txt',CoeffCldL,'\t');
CoeffAETL(isnan(CoeffAETL))=-9999;dlmwrite('CoeffNT1AETL.txt',CoeffAETL,'\t');
CoeffTmpL(isnan(CoeffTmpL))=-9999;dlmwrite('CoeffNT1TmpL.txt',CoeffTmpL,'\t');
CoeffCldG(isnan(CoeffCldG))=-9999;dlmwrite('CoeffNT1CldG.txt',CoeffCldG,'\t');
CoeffAETG(isnan(CoeffAETG))=-9999;dlmwrite('CoeffNT1AETG.txt',CoeffAETG,'\t');
CoeffTmpG(isnan(CoeffTmpG))=-9999;dlmwrite('CoeffNT1TmpG.txt',CoeffTmpG,'\t');
CoeffPercT1(isnan(CoeffPercT1))=-9999;dlmwrite('CoeffPercT1.txt',CoeffPercT1,'\t');
CoeffPercTmp(isnan(CoeffPercTmp))=-9999;dlmwrite('CoeffPercTmp.txt',CoeffPercTmp,'\t');
CoeffPercAET(isnan(CoeffPercAET))=-9999;dlmwrite('CoeffPercAET.txt',CoeffPercAET,'\t');
CoeffPercCld(isnan(CoeffPercCld))=-9999;dlmwrite('CoeffPercCld.txt',CoeffPercCld,'\t');

sigMonth(isnan(sigMonth))=-9999;dlmwrite('sigMonth.txt',sigMonth,'\t');

cIcld(isnan(cIcld))=-9999;dlmwrite('cIcld.txt',cIcld,'\t');
cIaet(isnan(cIaet))=-9999;dlmwrite('cIaet.txt',cIaet,'\t');
cItmp(isnan(cItmp))=-9999;dlmwrite('cItmp.txt',cItmp,'\t');
cIt1(isnan(cIt1))=-9999;dlmwrite('cIt1.txt',cIt1,'\t');

cI_coefcld(isnan(cI_coefcld))=-9999;dlmwrite('cI_coefcld.txt',cI_coefcld,'\t');
cI_coefaet(isnan(cI_coefaet))=-9999;dlmwrite('cI_coefaet.txt',cI_coefaet,'\t');
cI_coeftmp(isnan(cI_coeftmp))=-9999;dlmwrite('cI_coeftmp.txt',cI_coeftmp,'\t');
cI_coeft1(isnan(cI_coeft1))=-9999;dlmwrite('cI_coeft1.txt',cI_coeft1,'\t');

cIcld_S(isnan(cIcld_S))=-9999;dlmwrite('cIcld_S.txt',cIcld_S,'\t');
cIaet_S(isnan(cIaet_S))=-9999;dlmwrite('cIaet_S.txt',cIaet_S,'\t');
cItmp_S(isnan(cItmp_S))=-9999;dlmwrite('cItmp_S.txt',cItmp_S,'\t');
cIt1_S(isnan(cIt1_S))=-9999;dlmwrite('cIt1_S.txt',cIt1_S,'\t');

month_yr_stderror(isnan(month_yr_stderror))=-9999;dlmwrite('month_yr_stderror.txt',month_yr_stderror,'\t');