% Matlab-code for the paper entitled "Decreasing rainfall frequency
% contriutes to earlier leaf onset in northern ecosystems" 
% Authors: Jian Wang et al.
% Goal: determine climatic sensitivities using ridge regression
% T, P, I, W: temperature, precipitation amount, cloudiness, precipitation
% frequency
clc,clear
warning('off');
maindir1 = 'D:\data\LOD_NDVI3g\';
maindir2 = 'D:\data\preseason_vars\';
outdir = 'D:\result\sensitivity_NDVI3g\';
year = 34;
T_TOTAL = zeros(year, 720, 4320,'double');
P_TOTAL = zeros(year, 720, 4320,'double');
I_TOTAL = zeros(year, 720, 4320,'double');
W_TOTAL = zeros(year, 720, 4320,'double');
LOD_TOTAL = zeros(year, 720, 4320,'double');
Sen_T = zeros(720, 4320,'double');
Sen_P = zeros(720, 4320,'double');
Sen_I = zeros(720, 4320,'double');
Sen_W = zeros(720, 4320,'double');
P_value = zeros(720, 4320,'double');

for j = 1 : year
year_str = num2str(j+1981);
name_T = ['preseason_tmp_',year_str,'.tif'];
file_name_T = fullfile(maindir,name_T);
[T,info2] = geotiffread(file_name_T);
T_TOTAL(j,:,:) = T;
name_P = ['preseason_pre_',year_str,'.tif'];
file_name_P = fullfile(maindir,name_P);
[P,info3] = geotiffread(file_name_P);
P_TOTAL(j,:,:) = P;
name_I = ['preseason_cld_',year_str,'.tif'];
file_name_I = fullfile(maindir,name_I);
[I,info4] = geotiffread(file_name_I);
I_TOTAL(j,:,:) = I;
name_W = ['preseason_wet_',year_str,'.tif'];
file_name_W = fullfile(maindir,name_W);
[W,info5] = geotiffread(file_name_W);
W_TOTAL(j,:,:) = W;
name_LOD = ['LOD',year_str,'.tif'];
file_name_LOD = fullfile(maindir2,name_LOD);
[LOD,info6] = geotiffread(file_name_LOD);
LOD_TOTAL(j,:,:) = LOD;
end

for column = 1 : 4320
    for line = 1 : 720
        LOD_site = LOD_TOTAL(:,line,column);
        I_site = I_TOTAL(:,line,column);
        if length(unique(I_site)) < 17
            continue
        end
        T_site = T_TOTAL(:,line,column);
        P_site = P_TOTAL(:,line,column);
        W_site = W_TOTAL(:,line,column);
        
        LOD_s_anl = LOD_site - mean(LOD_site);
        LOD_s_anl_nor = (LOD_s_anl-min(LOD_s_anl))/(max(LOD_s_anl)-min(LOD_s_anl));
        LOD_s_a = LOD_s_anl_nor;
        T_s_anl = T_site - mean(T_site);
        T_s_anl_nor = (T_s_anl-min(T_s_anl))/(max(T_s_anl)-min(T_s_anl));
        T_s_a = T_s_anl_nor;
        P_s_anl = P_site - mean(P_site);
        P_s_anl_nor = (P_s_anl-min(P_s_anl))/(max(P_s_anl)-min(P_s_anl));
        P_s_a = P_s_anl_nor;
        I_s_anl = I_site - mean(I_site);
        I_s_anl_nor = (I_s_anl-min(I_s_anl))/(max(I_s_anl)-min(I_s_anl));
        I_s_a = I_s_anl_nor;
        W_s_anl = W_site - mean(W_site);
        W_s_anl_nor = (W_s_anl-min(W_s_anl))/(max(W_s_anl)-min(W_s_anl));
        W_s_a = W_s_anl_nor;
        
        x_a = [T_s_a,P_s_a,I_s_a,W_s_a];
        y_a = LOD_s_a;
        B_a = ridge(y_a,x_a,1,0);
        y_a_e = x_a*B_a(2:end)+B_a(1);
        [R_as,P_as] = corrcoef(y_a,y_a_e);
        Sen_T(line,column) = B_a(2);
        Sen_P(line,column) = B_a(3);
        Sen_I(line,column) = B_a(4);
        Sen_W(line,column) = B_a(5);
        P_value(line,column) = P_as(1,2);
    end
end
filename1 = fullfile(outdir,'sen_T.tif');
filename2 = fullfile(outdir,'sen_P.tif');
filename3 = fullfile(outdir,'sen_I.tif');
filename4 = fullfile(outdir,'sen_W.tif');
filename5 = fullfile(outdir,'P_value.tif');
geotiffwrite(filename1,Sen_T,info2);
geotiffwrite(filename2,Sen_P,info2);
geotiffwrite(filename3,Sen_I,info2);
geotiffwrite(filename4,Sen_W,info2);
geotiffwrite(filename5,P_value,info2);