%%These codes are provided to support paper 'Deceleration of China¡¯s human water use and its socioeconomic drivers' submitted to PNAS
%% developed by Yan Bo
%%This function script is supposed to run in Matlab ver. R2016b
%%Compatibility to other version of Matlab is not tested.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Part 1  piecewise regression revised based on Wang Xuhui 2012
%piecewise regression
%This function script is supposed to run in Matlab ver. 2010b
%Compatibility to other version of Matlab is not tested.
%PLEASE DO NOT DISTRIBUTE THE PROGRAM WITHOUT AUHTORS' PRIOR PERMISSION.
%Copyright Wang Xuhui 2012
%@Aug 2012, adding test for non-time series data
%INPUT PARAMETERS:
%st: earliest potential th
%ed: latest potential th
%OUTPUT:
%th: turning point
%p: p-stat for th
%mb,xx: data used for reconstruct the regression line outside the function
%p1: p-stat for the 1st part
%p2: p-stat for the 2nd part
%rsq: r-square of the model
%p_rsq: p-stat of the model fit

function [th,th2,p,mb,rsq] = pieceregWangForSun(y,x1,st,ed)

    if nargin<4
        error('not enough input arguments')
    end
    
    mres = 99999*10^5;
    myres = 99999*10^5;
    mpt = 0;
    yrlen = length(x1);
    mp = 1;
    
    for bpt = st:1:ed
        x2 = zeros(yrlen,1);
        aa = x1-bpt;
        x2(x1>bpt) = aa(x1>bpt);
        x2(x1<=bpt) = 0;
        
        for yr=1:yrlen
            if (yr>bpt)
                x2(yr) = yr-bpt;
            end
        end
        
        for bpt2 = bpt+1:1:ed
            x3 = zeros(yrlen,1);
            bb = x1-bpt2;
            x3(x1>bpt2) = bb(x1>bpt2);
            x3(x1<=bpt2) = 0;

            for yr=1:yrlen
                if (yr>bpt2)
                    x3(yr) = yr-bpt2;
                end
            end
            
            mystats = regstats(y,[x1,x2,x3],'linear',{'yhat','r','tstat'});
            yres = sum((mystats.r.*1000).^2)./1000000;
            if (yres < myres)
            x = [ones(yrlen,1),x1,x2,x3];
            [b,~,~,~,stats] = regress(y,x);
            myres = yres;
            mres = stats(4); 
            mpt = bpt; mpt2 = bpt2;
            mp = mystats.tstat.pval(3); %p value
            mb = b; %coeffs of reg
            xx = x;
            rsq = stats(1);
            p_rsq = stats(3);
            end
        end
    end
    
    th = mpt;
    th2 = mpt2;
    p = mp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Part 2
%% Data Process: Drivers of water use trends-LMDI results at three periods

[num,txt] = xlsread('LMDI.xlsx','Data','A3:BT16760');
% irrigation
y = num(:,3:7);
x1 = num(:,8);
x2 = num(:,9:13)./num(:,8);
x3 = num(:,15:19);

y(y==0)=0.00000001;
x1(x1==0)=0.00000001;
x2(x2==0)=0.00000001;
x3(x3==0)=0.00000001;

HWU = num(:,2);

ln_y = log(y);
ln_x1 = log(x1);
ln_x2 = log(x2);
ln_x3 = log(x3);

city = unique(txt(:,2));
year = 1956:1:2013;
year = year';
%a = [10:34;34:58];

S_HWU = [];
A = [];
B = [];
C = [];
D = [];
E=[];

CON_IRR = zeros(50,3);

for i = 1:length(city)
    row_index = find(strcmp(txt(:,2),city(i)));
    
    for k = 1:3
        if k==1
            a = [1 11];
        elseif k==2
            a = [11 28];
        elseif k==3
            a = [28 49];
        end

       
        S_IRR((i-1)*3+k,1) = HWU(row_index(a(2)))-HWU(row_index(a(1)));
        
        for j = 1:5
	A(j,1) = y(row_index(a(2)),j)-y(row_index(a(1)),j);
	B(j,1) = ln_y(row_index(a(2)),j)-ln_y(row_index(a(1)),j);
	C(j,1) = ln_x1(row_index(a(2)))-ln_x1(row_index(a(1)));
	D(j,1) = ln_x2(row_index(a(2)),j)-ln_x2(row_index(a(1)),j);
	E(j,1) = ln_x3(row_index(a(2)),j)-ln_x3(row_index(a(1)),j);
        end

        CON_IRR((i-1)*3+k,1) = nansum(A.*C./B,1);
        CON_IRR((i-1)*3+k,2) = nansum(A.*D./B,1);
        CON_IRR((i-1)*3+k,3) = nansum(A.*E./B,1);
    end
end



%% IND
y = num(:,21:31);
x1 = num(:,32);
x2 = num(:,33:43)./num(:,32);
x3 = num(:,45:55);

y(y==0)=0.00000001;
x1(x1==0)=0.00000001;
x2(x2==0)=0.00000001;
x3(x3==0)=0.00000001;

HWU = num(:,20);

ln_y = log(y);
ln_x1 = log(x1);
ln_x2 = log(x2);
ln_x3 = log(x3);


city = unique(txt(:,2));
year = 1956:1:2013;
year = year';


for i = 1:length(city)
    row_index = find(strcmp(txt(:,2),city(i)));
    
    for k = 1:3
        if k==1
            a = [1 11];
        elseif k==2
            a = [11 28];
        elseif k==3
            a = [28 49];
        end

        S_IND((i-1)*3+k,1) = HWU(row_index(a(2)))-HWU(row_index(a(1)));
        for j = 1:11
	A(j,1) = y(row_index(a(2)),j)-y(row_index(a(1)),j);
	B(j,1) = ln_y(row_index(a(2)),j)-ln_y(row_index(a(1)),j);
	C(j,1) = ln_x1(row_index(a(2)))-ln_x1(row_index(a(1)));
	D(j,1) = ln_x2(row_index(a(2)),j)-ln_x2(row_index(a(1)),j);
	E(j,1) = ln_x3(row_index(a(2)),j)-ln_x3(row_index(a(1)),j);
        end
        CON_IND((i-1)*3+k,1) = nansum(A.*C./B,1);
        CON_IND((i-1)*3+k,2) = nansum(A.*D./B,1);
        CON_IND((i-1)*3+k,3) = nansum(A.*E./B,1);
    end
end

%% Urban sector
y1 = num(:,57);
x11 = num(:,58);
x12 = num(:,59)*365;

y2 = num(:,60);
x21 = num(:,61);
x22 = num(:,62);

y1(y1==0)=0.00000001;
x11(x11==0)=0.00000001;
x12(x12==0)=0.00000001;

y2(y2==0)=0.00000001;
x21(x21==0)=0.00000001;
x22(x22==0)=0.00000001;

HWU = num(:,56);

ln_y1 = log(y1);
ln_x11 = log(x11);
ln_x12 = log(x12);

ln_y2 = log(y2);
ln_x21 = log(x21);
ln_x22 = log(x22);


city = unique(txt(:,2));
year = 1956:1:2013;
year = year';
%a = [10:34;35:58];


for i = 1:length(city)
    row_index = find(strcmp(txt(:,2),city(i)));
    
     for k = 1:3
        if k==1
            a = [1 11];
        elseif k==2
            a = [11 28];
        elseif k==3
            a = [28 49];
        end

        S_URB((i-1)*3+k,1) = HWU(row_index(a(2)))-HWU(row_index(a(1)));
        

      A = y1(row_index(a(2)))-y1(row_index(a(1)));
      B = ln_y1(row_index(a(2)))-ln_y1(row_index(a(1)));
      C = ln_x11(row_index(a(2)))-ln_x11(row_index(a(1)));
      D = ln_x12(row_index(a(2)))-ln_x12(row_index(a(1)));

       E = y2(row_index(a(2)))-y2(row_index(a(1)));
       F = ln_y2(row_index(a(2)))-ln_y2(row_index(a(1)));
       G = ln_x21(row_index(a(2)))-ln_x21(row_index(a(1)));
       H = ln_x22(row_index(a(2)))-ln_x22(row_index(a(1)));
        
       CON_URB((i-1)*3+k,1) = nansum(A.*C./B,1);
       CON_URB((i-1)*3+k,2) = nansum(A.*D./B,1);
       CON_URB((i-1)*3+k,3) = nansum(E.*G./F,1);
       CON_URB((i-1)*3+k,4) = nansum(E.*H./F,1);
    end
end

%% rural sector
y1 = num(:,64);
x11 = num(:,65);
x12 = num(:,66)*365;

y2 = num(:,67);
x21 = num(:,68);
x22 = num(:,69)*365;

y1(y1==0)=0.00000001;
x11(x11==0)=0.00000001;
x12(x12==0)=0.00000001;

y2(y2==0)=0.00000001;
x21(x21==0)=0.00000001;
x22(x22==0)=0.00000001;

HWU = num(:,63);

ln_y1 = log(y1);
ln_x11 = log(x11);
ln_x12 = log(x12);

ln_y2 = log(y2);
ln_x21 = log(x21);
ln_x22 = log(x22);

city = unique(txt(:,2));
year = 1956:1:2013;
year = year';
%a = [10:34;34:58];

CON_RUR = zeros(50,3);
for i = 1:length(city)
    row_index = find(strcmp(txt(:,2),city(i)));
    
   for k = 1:3
        if k==1
            a = [1 11];
        elseif k==2
            a = [11 28];
        elseif k==3
            a = [28 49];
        end

        S_RUR((i-1)*3+k,1) = HWU(row_index(a(2)))-HWU(row_index(a(1)));
        A = y1(row_index(a(2)))-y1(row_index(a(1)));
        B = ln_y1(row_index(a(2)))-ln_y1(row_index(a(1)));
        C =ln_x11(row_index(a(2)))-ln_x11(row_index(a(1)));
        D = ln_x12(row_index(a(2)))-ln_x12(row_index(a(1)));

        E = y2(row_index(a(2)))-y2(row_index(a(1)));
        F = ln_y2(row_index(a(2)))-ln_y2(row_index(a(1)));
        G = ln_x21(row_index(a(2)))-ln_x21(row_index(a(1)));
        H = ln_x22(row_index(a(2)))-ln_x22(row_index(a(1)));
        
        CON_RUR((i-1)*3+k,1) = nansum(A.*C./B,1);
        CON_RUR((i-1)*3+k,2) = nansum(A.*D./B,1);
        CON_RUR((i-1)*3+k,3) = nansum(E.*G./F,1);
        CON_RUR((i-1)*3+k,4) = nansum(E.*H./F,1);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Part 3
%% structural process-based model

%% IRR
clear all
clc
[num,txt] = xlsread('E:\WU paper\Contribution analysis.xlsx','IRR-reg');
txt(1,:) = [];
Pro_ID = unique(txt(:,1));
modelfun = @(b,x)(x(:,1).*(b(1).*x(:,2)+b(2)));
for i = 1:length(Pro_ID)
    row_index = find(strcmp(Pro_ID(i),txt(:,1))==1&num(:,1)>=1975&num(:,1)<=2010);
    beta0 = [0.01 0.01];
    [x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(modelfun,beta0,[num(row_index,6),num(row_index,9)],num(row_index,11),[-inf -inf],[0 inf]);
    beta(i,:) = x;
end


%% IND
[num,txt] = xlsread('E:\WU paper\Contribution analysis.xlsx','IND-reg','A1:F1569');
txt(1,:) = [];
Pro_ID = unique(txt(:,1));
modelfun = @(b,x)b(1).*((1-(1-x(:,1)).*x(:,2)))./(1-b(2).*((1-(1-x(:,1)).*x(:,2))));

for i = 1:length(Pro_ID)
    row_index = find(strcmp(Pro_ID(i),txt(:,1))==1&num(:,1)>=1975&num(:,1)<=2013);
    b = 0.1;
    a = 0.1;    
    beta0 = [a b];
    [x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(modelfun,beta0,[num(row_index,3),num(row_index,2)],num(row_index,5),[0 0],[inf 1]);
    beta(i,:) = x;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Part 4
%% Figure 1

clear all
clc
[num,txt] = xlsread('chaningpoint.xlsx','Data');
[num2,txt2] = xlsread('chaningpoint.xlsx','Area');
txt2(1,:) = [];
x = 1965:1:2013;
x = x';

%%random resampling
for i = 1:300
    disp(i)
    ID = randperm(341,ceil(341*0.90));
    ID = sort(ID);
   
    for j = 1:length(ID)
        index = find(strcmp(txt(1,ID(j)+1),txt2(:,2))==1);
        Area_sample(j) = num2(index,4);
    end
    WU_sample = sum(num(:,ID+1),2);
    WU_sample_mm = WU_sample./sum(Area_sample).*1000000;

    [th,th2,p,mb,rsq]  = CP_REV_BY(WU_sample_mm,x,1965,2013);

    year(i,1) = th;
    year(i,2) = th2;
    slope(i,1) = mb(2);
    slope(i,2) = mb(2)+mb(3);
    slope(i,3) = mb(2)+mb(3)+mb(4);
   
    intercept(i,1) = WU_sample_mm(th-1965+1)-slope(i,1)*th;
    intercept(i,2) = WU_sample_mm(th-1965+1)-slope(i,2)*th;%num(th2-1965+1,end)-s2*th2;
    intercept(i,3) = slope(i,2)*th2+intercept(i,2)-slope(i,3)*th2;%WU_sample_mm(th2-1965+1)-slope(i,3)*th2;WU_sample_mm(th-1965+1)-slope(i,2)*th2-slope(i,3)*th2
   
   rsq_final(i,1) = rsq;
   rsq_final_ad(i,1) = 1-(1-rsq)*(49-1)/(49-4);
end

%% Figure 1a
figure('position',[100 100 700 350])
WU_mm = num(:,end)./sum(num2(:,4))*1000000;

%frequency distribution of changing point 1
[f,xi] = ksdensity(year(:,1));
index1 = find(cumsum(f)>=sum(f)*0.25&cumsum(f)<=sum(f)*0.75);
CT1 = cbrewer('seq','Blues',length(index1));
[B,I] = sort(f(index1));
for i = 1:length(index1)
    %plot([xi(I(i)) xi(I(i))],[21 79],'-','Color',CT1(i,:),'LineWidth',2)
    h1 = patch([xi(index1((I(i))))-0.29 xi(index1((I(i))))+0.29 xi(index1((I(i))))+0.29 xi(index1((I(i))))-0.29],[20.5 20.5 74.5 74.5],'-','FaceColor',CT1(i,:),'FaceAlpha',0.6,'EdgeColor','none');
   
    hold on
end

%frequency distribution of changing point 1
[f,xi] = ksdensity(year(:,2));
index2 = find(cumsum(f)>=sum(f)*0.25&cumsum(f)<=sum(f)*0.75);
CT2 = cbrewer('seq','Oranges',length(index2));
[B,I] = sort(f(index2));
for i = 1:length(index2)
    %plot([xi(I(i)) xi(I(i))],[21 79],'-','Color',CT1(i,:),'LineWidth',2)
    h1 = patch([xi(index2((I(i))))-0.2 xi(index2((I(i))))+0.2 xi(index2((I(i))))+0.2 xi(index2((I(i))))-0.2],[20.5 20.5 74.5 74.5],'-','FaceColor',CT2(i,:),'FaceAlpha',0.6,'EdgeColor','none');
    hold on
end
hold on

%plot changing point
plot([1975 1975],[20 75],'--','Color','k','LineWidth',1.3)
plot([1992 1992],[20 75],'--','Color','k','LineWidth',1.3)

% grey lines
for i = 1:300
    x1 = [1965 year(i,1)];
    x2 = [year(i,1) year(i,2)];
    x3 = [year(i,2) 2013];
    plot(x1,intercept(i,1)+slope(i,1)*x1,'-','Color',[230 230 230]/255,'LineWidth',1.3);
    plot(x2,intercept(i,2)+slope(i,2)*x2,'-','Color',[230 230 230]/255,'LineWidth',1.3);
    plot(x3,intercept(i,3)+slope(i,3)*x3,'-','Color',[230 230 230]/255,'LineWidth',1.3);
    hold on
end


h = plot(1965:1:2013',WU_mm,'Marker','o','MarkerSize',4)
% h = plot(1965:1:2013',num(:,end),'Marker','o','MarkerSize',4)
h.MarkerFaceColor = 'k'
h.Color = 'k'
h.LineWidth = 1.3
[th,th2,p,mb,rsq]  = CP_REV_BY(WU_mm,x,1965,2013);
s1 = mb(2);
s2 = mb(2)+mb(3);
s3 = mb(2)+mb(3)+mb(4); 
inte1 = WU_mm(th-1965+1)-s1*th;
inte2 = WU_mm(th-1965+1)-s2*th;%num(th2-1965+1,end)-s2*th2;
inte3 = s2*th2+inte2-s3*th2;
hold on
x1 = [1965 th];
x2 = [th th2];
x3 = [th2 2013];
plot(x1,inte1+s1*x1,'-','Color','r','LineWidth',1.3)
plot(x2,inte2+s2*x2,'-','Color','r','LineWidth',1.3)
plot(x3,inte3+s3*x3,'-','Color','r','LineWidth',1.3)

xlabel({'Year'})
ylabel({'Water use (mm yr^{-1})'})
set(gca,'FontSize',12,'linewidth',1.2,'FontName','Arial');
xlim([1964 2014])
xticks([1965:8:2013])
ylim([26 75])
yticks([30:10:70])
box on

file_name = ['FIG1_CPv2.tif'];
print(gcf, '-r600', '-dtiff',file_name);


%%Figure 1b
figure('position',[100 100 400 300])
data = [s1,s2,s3;median(slope(:,1)),median(slope(:,2)),median(slope(:,3))];
data(3,:) = data(2,:);
data(2,:) = nan;
h = bar(data');

color = [0 0 0;255 255 255;230 230 230]/255
for i = 1:3
    h(i).EdgeColor = 'k';
    h(i).LineWidth = 1.0;
    h(i).FaceColor = color(i,:);
    h(i).BarWidth = 1.6;
end
hold on
data2 = [prctile(slope(:,1),25),prctile(slope(:,2),25),prctile(slope(:,3),25);prctile(slope(:,1),75),prctile(slope(:,2),75),prctile(slope(:,3),75)];
for i = 1:3
    t = errorbar(0.2+i,data(3,i),data2(1,i)-data(3,i),data2(2,i)-data(3,i))
    t.Color = 'k';
    t.LineWidth = 1;  
end

xticklabels({'P1','P2','P3'})
ylabel({'Water use trend (mm yr^{-2})'})
ylim([0.2 1.3])
xlim([0.5 3.5])
set(gca,'FontSize',14,'linewidth',1.2,'FontName','Arial');

file_name = ['FIG12_CP_slope.tif'];
print(gcf, '-r600', '-dtiff',file_name);


%% Fugure 1c
[num,txt] = xlsread('E:\WU paper\Trend.xls','slown_down_pattern');
txt(1,:) = [];
P = unique(txt(:,4));
C = [146 206 80;255 192 0]/255
figure('position',[100 100 400 328])
for i = 1:length(P)
    row_index = find(strcmp(txt(:,4),P(i))==1);
    x_min = min(num(row_index,11));
    x_max = max(num(row_index,11));
    x_ran = x_min:0.1:x_max;
    for j = 1:length(x_ran)
        row_index2 = find(num(row_index,11)<=x_ran(j));
        y(j) = sum(num(row_index(row_index2),9))./sum(num(row_index,9))*100;
    end
    h(i) = plot(x_ran,y,'LineWidth',2.3)
    x_ran = [];
    y = [];    
    hold on
   
end

legend([h(1),h(2)],{'P2 v.s. P1','P3 v.s. P2'},'box','off','location','northwest')
hold on
set(gca,'FontSize',14,'linewidth',1.2,'FontName','Arial');

plot([0 0],[0 100],'--','Color','k','LineWidth',1.3)
xlabel({'Difference in WU trend (mm yr^{-2})'})
ylabel({'Cumulative proportion (%)'})
xlim([-10 10])

file_name = ['FIG1.3_CP_slopev2.tif'];
print(gcf, '-r600', '-dtiff',file_name);

%% Figure 3
clear all
[num,txt] = xlsread('LMDI.xlsx','results_3p','B1:Q4');
CT = [9 81 157;33 114 180;66 146 197;64 171 93;114 197 115;161 217 156;170 13 24;239 58 47;251 105 80;253 146 112;84 40 136;107 82 161;127 126 183;156 155 199]/255;
do = num(:,2:15);

figure('position',[100 100 1000 800]);
axes('position',[0.17 0.15 0.8 0.8]);

s=0;
    l=1;
    flag = 1;

    ggg = [10,17,21];
for i = 1:3

    b=0;
    c = num(i,1);
    
    index1 = find(do(i,:)>=0);
    [B1,I1] = sort(do(i,index1));
    index2 = find(do(i,:)<0);
    [B2,I2] = sort(do(i,index2),'descend');
    do2 = [do(i,index1(I1)),do(i,index2(I2))];
    a = [s+l*(flag-1) s+l*flag s+l*flag s+l*(flag-1)];
    if i==1
    patch(a,[b b c c],[150 150 150]/255,'EdgeColor','none');%1965
    end
    flag = flag+1;
    b = b+num(i,1);
    for j = 1:14
        c = c+do2(j);
        a = [s+l*(flag-1) s+l*flag s+l*flag s+l*(flag-1)];
        patch(a,[b b c c],CT(find(do2(j)==do(i,:)),:),'EdgeColor','none');%1965
        aaaa(i,j) = find(do2(j)==do(i,:));
        flag = flag+1;
        if i==1|i==2
            text(0.5+j+15*(i-1),max(b,c)+15,txt(1,find(do2(j)==do(i,:))+1),'rotation',90,'Color',CT(find(do2(j)==do(i,:)),:),'FontName','Arial','FontSize',10);
            con = num(i,find(do2(j)==do(i,:))+1)/num(i,1)*100*100;
            con = round(con);
            con = con/100;            
            text(0.5+j+15*(i-1),min(b,c)-70,[num2str(con),'%'],'rotation',90,'Color',CT(find(do2(j)==do(i,:)),:),'FontName','Arial','FontSize',10);
        elseif i==3&j<9
            text(0.5+j+15*(i-1),max(b,c)+15,txt(1,find(do2(j)==do(i,:))+1),'rotation',90,'Color',CT(find(do2(j)==do(i,:)),:),'FontName','Arial','FontSize',10);

            con = num(i,find(do2(j)==do(i,:))+1)/num(i,1)*100*100;
            con = round(con);
            con = con/100;        
            text(0.5+j+15*(i-1),min(b,c)-70,[num2str(con),'%'],'rotation',90,'Color',CT(find(do2(j)==do(i,:)),:),'FontName','Arial','FontSize',10);
        else
            text(0.5+j+15*(i-1),min(b,c)-13,txt(1,find(do2(j)==do(i,:))+1),'rotation',270,'Color',CT(find(do2(j)==do(i,:)),:),'FontName','Arial','FontSize',10);

            con = num(i,find(do2(j)==do(i,:))+1)/num(i,1)*100*100;
            con = round(con);
            con = con/100;    
            text(0.5+j+15*(i-1),max(b,c)+10,[num2str(con),'%'],'rotation',90,'Color',CT(find(do2(j)==do(i,:)),:),'FontName','Arial','FontSize',10);
        end
        b = b+do2(j);
    end
    a = [s+l*(flag-1) s+l*flag s+l*flag s+l*(flag-1)];
    patch(a,[0 0 num(i,end) num(i,end)],[150 150 150]/255,'EdgeColor','none');%1965 
end

   ylim([200 1005]);
    xlim([0 46])
xticks([0:1:46])
xticklabels({ })

ylabel({'Water use (km^{3} yr^{-1})'});
box on
set(gca,'FontName','Arial','FontSize',15,'linewidth',1.2);

labels = [1965,1975,1992,2013];
text(-0.8,200-30,num2str(labels(1)),'FontName','Arial','FontSize',15)
text(15-0.8,200-30,num2str(labels(2)),'FontName','Arial','FontSize',15)
text(30-0.8,200-30,num2str(labels(3)),'FontName','Arial','FontSize',15)
text(45-0.8,200-30,num2str(labels(4)),'FontName','Arial','FontSize',15)

file_name = ['3.LMDIv3.tif'];
print(gcf, '-r600', '-dtiff',file_name);

%%Figure 4
%% figure 4a&b IRR
color = [47 137 190;103 194 161;127 126 184;180 180 180;]/255
[num,txt] = xlsread('Data.xlsx','Contribution analysis','B3:H33');
[num2,txt2] = xlsread('Data.xlsx','Contribution analysis','L3:R33');

for i = 1:2
    if i==1
        X = num(:,3:end);

    else
        X = num2(:,3:end);

    end
    Xneg = X;
    Xneg(Xneg>=0)=0;
    Xpos = X;
    Xpos(Xpos<0)=0;
    
    figure('position',[100 100 1000 300])
    a = axes('position',[0.15 0.3 0.8 0.66])
    h = bar([1:1:31],Xneg,'stacked');
    xlim([0 32])
    for c = 1:4
        h(c).FaceColor = color(c,:);
        h(c).EdgeColor = 'none';
    end
    hold on
    h = bar([1:1:31],Xpos,'stacked');

    for c = 1:4
        h(c).FaceColor = color(c,:);
        h(c).EdgeColor = 'none';
    end
    xticks([1:1:31])
    if i==1
        xticklabels(txt);
    else
        xticklabels(txt2);
    end
    set(gca,'XTickLabelRotation',35);
    hold on
    if i==1
    a = scatter([1:1:31],num(:,1),10,'MarkerFaceColor','k','MarkerEdgeColor','k')
    b = scatter([1:1:31],num(:,2),10,'MarkerFaceColor',[211 63 77]/255,'MarkerEdgeColor',[211 63 77]/255)
    else
    a = scatter([1:1:31],num2(:,1),10,'MarkerFaceColor','k','MarkerEdgeColor','k')
    b = scatter([1:1:31],num2(:,2),10,'MarkerFaceColor',[211 63 77]/255,'MarkerEdgeColor',[211 63 77]/255)
    end
    ylabel([{'Trend in WUI'},{'(mm yr^{-2})'}])
    
    
    set(gca,'FontName','Arial','FontSize',13,'linewidth',1.2);
    
    legend([a,b,h(1:4)],{'OBS','SIM','WCI','PIR','AIR','OF'},'box','off','Orientation','Horizontal','location','southeast');
    
    file_name = ['4.IRR_Drivers_P',num2str(i+1),'v3.tif'];
    print(gcf, '-r600', '-dtiff',file_name);
end

%% figure 4c&d IND
color = [47 137 190;103 194 161;180 180 180;]/255;
[num,txt] = xlsread('Data.xlsx','Contribution analysis','B39:G69');
[num2,txt2] = xlsread('Data.xlsx','Contribution analysis','J39:O69');

for i = 1:2
    if i==1
        X = num(:,3:end);

    else
        X = num2(:,3:end);

    end
    Xneg = X;
    Xneg(Xneg>=0)=0;
    Xpos = X;
    Xpos(Xpos<0)=0;
    
    figure('position',[100 100 1000 300])
    a = axes('position',[0.15 0.3 0.8 0.66])
    h = bar([1:1:31],Xneg,'stacked');
    xlim([0 32])
    for c = 1:3
        h(c).FaceColor = color(c,:);
        h(c).EdgeColor = 'none';
    end
    hold on
    h = bar([1:1:31],Xpos,'stacked');

    for c = 1:3
        h(c).FaceColor = color(c,:);
        h(c).EdgeColor = 'none';
    end
    xticks([1:1:31])
    if i==1
        xticklabels(txt);
    else
        xticklabels(txt2);
    end
    set(gca,'XTickLabelRotation',35);
    hold on
    if i==1
    a = scatter([1:1:31],num(:,1),10,'MarkerFaceColor','k','MarkerEdgeColor','k')
    b = scatter([1:1:31],num(:,2),10,'MarkerFaceColor',[211 63 77]/255,'MarkerEdgeColor',[211 63 77]/255)
    else
    a = scatter([1:1:31],num2(:,1),10,'MarkerFaceColor','k','MarkerEdgeColor','k')
    b = scatter([1:1:31],num2(:,2),10,'MarkerFaceColor',[211 63 77]/255,'MarkerEdgeColor',[211 63 77]/255)
    end
     ylabel([{'Trend in WUI'},{'(m^{3} Yuan^{-1} yr^{-2})'}])
    
    set(gca,'FontName','Arial','FontSize',13,'linewidth',1.2);
    
    legend([a,b,h(1:3)],{'OBS','SIM','REC','EVA','OF'},'box','off','Orientation','Horizontal','location','southeast');
    
    file_name = ['4.IND_Drivers_P',num2str(i+1),'v3.tif'];
    print(gcf, '-r600', '-dtiff',file_name);
end

%% END