%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This script is used to plot Figure 1 reported in:                                  
%                                                                              
% Daniel Nohrstedt, Maurizio Mazzoleni, Charles F. Parker, and Giuliano Di
% Baldassarre. Policy change after natural hazard
% events, Nature Communication
%                                                                              
% This script was prepared by Maurizio Mazzoleni                 
% Uppsala University, 2020                                                  
% mail to: maurizio.mazzoleni@geo.uu.se 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear 
close all

Country_code=1;             %%Plot country code: Use 1 to plot country code or 0 otherwise
Norm_PFAchange=0;           %%Plot normalized PFA changes if =1, otherwise plot results of figure 1
Hazard={'All'};


%%Define anlaysis period for the estimation of the normalized hazard indices
lower_limit=2007;
higher_limit=2011;
lower_limit_baseline=1980;
higher_limit_baseline=2011;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Run pre-processing scripts
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
run load_data
run processing_policy_data
run processing_hazard_data


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Run data analysis scripts
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
run merge_policy_hazard
run analysis_baseline30yr


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Plot Figure 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xlim_max=[9,9,9,2.2];
xlim_min=[0,0,0,0.7];

%%Color scatter plots divided by income class
color_L=[215 48 39]./255;
color_LM=[252 141 89]./255;
color_UM=[254 225 144]./255;
color_H=[69 117 180]./255;
color_NaN=[0 0 0]./255;


%%Label x-axis
x1='Fatalities (normalized)';
x2='Affected people (normalized)';
x3='Economic losses (normalized)';
x4='Number of events (normalized)';
xlabel_text={x1,x2,x3,x4};
    
%%Label y-axis
if Norm_PFAchange==0
    ylabel_text='Policy changes (aggregated PFAs)'; 
elseif Norm_PFAchange==1
    ylabel_text='Policy changes (normalized PFAs)'; 
end
        
   
vec_country_plot=analysis_country(:,1)';
    
%%Define image quality and name
ImageFontSize=14;
LabelFontSize=14;
scatter_size=50;
LineWidth_zero=1;
AxisFontName='Arial';
    
    
for p=1:4
    figure
    %%Define input for the scatter plot
    x = (Analysis_mean(:,p));
    y = (Analysis_mean(:,5));
    cls = Income_2015(:,2);
    %%Split the data in different income groups depending on the value in cls
    x_L = x(cls==1);
    y_L = y(cls==1);
    x_LM = x(cls==2);
    y_LM = y(cls==2);
    x_UM = x(cls==3);
    y_UM = y(cls==3);
    x_H = x(cls==4);
    y_H = y(cls==4);
    x_NaN = x(cls==5);
    y_NaN = y(cls==5);

    %%Parameter horintal line
    x_line1=0:12;
    y_line1=0;
    x_line2=1;
    y_line2=-5:5;

    plot(x_line1,y_line1*ones(size(x_line1)),'--k','LineWidth',LineWidth_zero)
    hold on
    plot(x_line2*ones(size(y_line2)),y_line2,'--k','LineWidth',LineWidth_zero)
    hold on
    scatter(x_L,y_L,scatter_size,color_L,'filled','MarkerEdgeColor','k')
    hold on
    scatter(x_LM,y_LM,scatter_size,color_LM,'filled','MarkerEdgeColor','k')
    hold on
    scatter(x_UM,y_UM,scatter_size,color_UM,'filled','MarkerEdgeColor','k')
    hold on
    scatter(x_H,y_H,scatter_size,color_H,'filled','MarkerEdgeColor','k')
    hold on
    scatter(x_NaN,y_NaN,scatter_size,color_NaN,'filled','MarkerEdgeColor','k')
    hold on
    if Country_code==1  %%Loop to assign country code to each point of the scatter plot
        for iji=1:size(Analysis_0711,2)
            for rr=1:size(vec_country_plot,2)
                if strcmp(analysis_country(iji,1),vec_country_plot{1,rr})==1
                    if ~isnan(x(iji,1)) || ~isnan(x(iji,1))
                        ta = text(x(iji,1),y(iji,1),analysis_country(iji,1));
                        ta.FontSize = 8;
                    end
                    hold on
                end
            end
            hold on
        end
    end
    xlim([xlim_min(p) xlim_max(p)])
    ylim([-1.5 1.6])
    grid on
    xlabel(xlabel_text{p},'FontSize',LabelFontSize)
    ylabel('Average PFA changes','FontSize',LabelFontSize)
    set(gca,'FontName',AxisFontName,'FontSize',ImageFontSize)
    set(gca,'children',get(gca,'children')) %Send to back the ensemble
    ytickformat('%.3g')
    xtickformat('%.2g')

end