clc
clear
close all

addpath(genpath('')); 

cd('');

%% Load data

% water_scarcity_05dgr - Water Scarcity Index (WCI & WSI combined)
% gdp_per_capita_05dgr - GDP (US$) per capita (IIASA 2007)
% hum_ftprint_05dgr - Human Footprint Index, the relative human influence in
%  terrestrial biomes (0 = least influenced parts, 100 = most influenced
%  parts of the biome)
%  (http://sedac.ciesin.columbia.edu/data/set/wildareas-v2-human-footprint-geographic)
% natural_haz_05dgr - global multihazard frequency and distribution
%  (http://sedac.ciesin.columbia.edu/data/set/ndh-multihazard-frequency-distribution)
% hdi_05dgr - Human Development Index
% (UNDP: http://hdr.undp.org/en/composite/HDI)

load data_lecture6.mat

%% plot the rasters

figure;
subplot(2,3,1)
imagesc(gdp_per_capita_05dgr);
colorbar
title('Economics - GDP/cap');

subplot(2,3,2)
imagesc(hum_ftprint_05dgr);
colorbar
title('Human Footprint')

subplot(2,3,3)
imagesc(natural_haz_05dgr);
colorbar
title('Natural Hazards')

subplot(2,3,4);
imagesc(water_scarcity_05dgr);
set(gca, 'CLim', [0, 1]);
colorbar
title('Water Stress Index - WSI')

subplot(2,3,5);
imagesc(hdi_05dgr);
set(gca, 'CLim', [0, 1]);
colorbar
title('Human Development Index - HDI')

%% read country data

% fao_id - compatible with the adm0_05dgr 
% Region id - geographical region
% polity - governance level - the higher the more open governance (The
%  Worldwide Governance Indicators, 2015 Update Aggregate Governance
%  Indicators 1996-2014 www.govindicators.org

raw_country = readtable('governance_index.xlsx'); %read xls raw data

governance_indicator(:,1) = raw_country.FAO_ID;
governance_indicator(:,2) = raw_country.RegionID;
governance_indicator(:,3) = raw_country.governanceEffectiveness;

%% Change NaN in governance data to regional average

governance_indicator_reg(:,1) = unique(governance_indicator(:,2));

temp = accumarray(governance_indicator(:,2),governance_indicator(:,3),[],@nanmean); 
governance_indicator_reg(:,2) = temp;
clear temp

% fill NaN values in governance_indicator with regional averages
for j = 1:size(governance_indicator,1)
    if isnan(governance_indicator(j,3))
        governance_indicator(j,3) = governance_indicator_reg(governance_indicator(j,2),2);
    end
end

clear governance_indicator_reg

%% Governance data to raster

load('clip_zones.mat','adm0_05dgr'); %load country raster

%create an empty matrix with the same row and col dimensions as adm0_05dgr
governance_05dgr = zeros(size(adm0_05dgr)); 

for i = 1:size(governance_indicator,1) %loop that goes through all the countries listed in governance_indicator
    cntry_logical = adm0_05dgr == governance_indicator(i,1); %returns 1 whenever adm0 has the country id on row i
    governance_05dgr(cntry_logical) = governance_indicator(i,3);%collect the country data in the temp raster
end

governance_05dgr(governance_05dgr == 0) = NaN;

% plot values
figure;
imagesc(governance_05dgr);
colorbar
title('Governance');

%% Transform GDP/cap to logarithmic scale
% observed to work better for this kind of analysis

gdp_per_capita_05dgr_log = log(gdp_per_capita_05dgr);

%% Vulnerability and Resilience Index for each grid cell

% Resilience index is calculated based on 3 societal indicators: 
% Governance, GDP/cap & HDI
% -> Create an empty 3-layered variable
res_parameters = zeros(size(adm0_05dgr,1),size(adm0_05dgr,2),3);

% A. collect the indicators:

% 1. governance
res_parameters(:,:,1) = governance_05dgr;

% 2. economics: GDP per capita
res_parameters(:,:,2) = gdp_per_capita_05dgr_log;

% 3. social capital: HDI
res_parameters(:,:,3) = hdi_05dgr;


% Vulnerability index is calculated based on 3 physical indicators: 
% Water Scarcity, Natural Hazards & Human Footprint
% -> Create an empty 3-layered variable
vuln_parameters = zeros(size(adm0_05dgr,1),size(adm0_05dgr,2),3);

% 1. human footprint
vuln_parameters(:,:,1) = hum_ftprint_05dgr;

% 2. natural hazards
vuln_parameters(:,:,2) = natural_haz_05dgr;

% 3. water scarcity
vuln_parameters(:,:,3) = water_scarcity_05dgr;

% B. Normalisation by Scaling Between 0 and 1 
% normalised(e_i) = (e_i - Emin)/(Emax-Emin)
% some datasets have really high or low values, and those cannot be
% included into normalisation, thus we use the 5th and 95th percentile 
% instead of min and max

res_parameters_norm = zeros(size(res_parameters));
vuln_parameters_norm = zeros(size(vuln_parameters));

for i = 1:3
    temp_res = res_parameters(:,:,i);
    temp_vuln = vuln_parameters(:,:,i);
    
    % normalise resilience index
    if max(max(temp_res)) > 1
        temp_5perc = prctile(temp_res(:),5);
        temp_95perc = prctile(temp_res(:),95);
        res_index = (temp_res - temp_5perc) ./ (temp_95perc-temp_5perc);
    else
        res_index = temp_res;
    end
    
    % normalise vulnerability index
    if max(max(temp_vuln)) > 1
        temp_5perc = prctile(temp_vuln(:),5);
        temp_95perc = prctile(temp_vuln(:),95);
        vuln_index = (temp_vuln - temp_5perc) ./ (temp_95perc-temp_5perc);
    else
        vuln_index = temp_vuln;
    end
    
    res_index(res_index > 1) = 1;
    res_index(res_index < 0) = 0;
    vuln_index(vuln_index > 1) = 1;
    vuln_index(vuln_index < 0) = 0;
    
    
    res_parameters_norm(:,:,i) = res_index;
    vuln_parameters_norm(:,:,i) = vuln_index;

end

% calculate vulnerability and resilience indices
vuln_index = nanmean(vuln_parameters_norm,3);
res_index = nanmean(res_parameters_norm,3);

%% plot the parameters and the vulnerability index

% parameters
title_res_parameters = {'Governance','Economy (GDP/cap)','Social capital (HDI)'}; % NOTE: curly brackets needed for string type array
title_vuln_parameters = {'Human Footprint','Natural Hazards','Water Scarcity'}; % NOTE: curly brackets needed for string type array

% Check that the scales are ok:
figure
for i = 1:3
    subplot(2,3,i)
    imagesc(res_parameters_norm(:,:,i))
    title(title_res_parameters{i})
    colorbar
    subplot(2,3,i+3)
    imagesc(vuln_parameters_norm(:,:,i))
    title(title_vuln_parameters{i})
    colorbar
end


% index
figure;
    subplot(1,2,1)
    imagesc(res_index);
    colorbar
    title('resilience index');
    subplot(1,2,2)
    imagesc(vuln_index);
    colorbar
    title('vulnerability index');

    
%% Calculate the individual vulnerability and resilience parameters for one basin

% for one basin, e.g. Ganges-Brahmaputra (id = 7)

load('clip_zones.mat', 'basins_05dgr');

% resilience parameters need to be weighted with population
% vulnerability parameters need to be weighted with surface area

load('total_pop.mat', 'total_pop_sel');

load('support_05dgr.mat','area_05dgr');

basin_logical = basins_05dgr == 7;

% For resilience parameters calculate the sum of population weighted parameters
% inside ganges basin and divide with total population of the basin
    for i = 1:3
        res_parameters_norm_weighted(:,:,i) = res_parameters_norm(:,:,i) .* total_pop_sel(:,:,5);
        temp_weighted = basin_logical .* res_parameters_norm_weighted(:,:,i);  
        temp_pop = basin_logical .* total_pop_sel(:,:,5);
        ganges_res(1,i) = nansum(nansum(temp_weighted))./nansum(nansum(temp_pop));
    end

% For vulnerability parameters calculate the sum of area weighted parameters
% inside ganges basin and divide with total area of the basin    
    for i = 1:3
        vuln_parameters_norm_weighted(:,:,i) = vuln_parameters_norm(:,:,i) .* area_05dgr;
        temp_weighted = basin_logical .* vuln_parameters_norm_weighted(:,:,i);    
        temp_area = basin_logical .* area_05dgr;
        ganges_vuln(1,i) = nansum(nansum(temp_weighted))./nansum(nansum(temp_area));
    end
clear temp*

ganges_res(1,4) = nanmean(ganges_res);    
ganges_vuln(1,4) = nanmean(ganges_vuln);    

%% EXAMPLE: Plot the results as a bar chart

% bar chart of different parameters
figure
title_vuln = [title_vuln_parameters,'total'];
barh(ganges_vuln, 0.5); % 0.5 is the relative bar width 
set(gca,'YTickLabel',title_vuln);

% stacked bar chart of vuln index
ganges_vuln_stacked(1,1:3) = zeros; %stacked bar chart only works with 2 or more stacked bars, so we need to fool matlab by creating an extra dataset
ganges_vuln_stacked(2,:) = ganges_vuln(1:3) / 3;

figure
barh(ganges_vuln_stacked, 0.5, 'stacked');
set(gca,'Ylim',[1.5 2.5], 'YTick', 2, 'YTickLabel', 'vulnerability index'); 
%if we set the y-axis limit like this, we don't see the extra dataset. In
%this case we also need to specify where the label should be by using YTick

% combine the two charts
new_label_vuln = [title_vuln(4) title_vuln(1:3)];

figure 
barh([1 2], ganges_vuln_stacked, 0.5, 'stacked');
hold on
h = barh(3:5, ganges_vuln(1:3), 0.5);
set(gca,'Ylim',[1.5 5.5], 'YTick', 2:8,'YTickLabel',new_label_vuln, 'Ticklength', [0,0]);

% if you want the parameters to have the same colors as in the stacked bar,
% it gets a bit tricky but can be done:
% the colors of the 5 parameters will be different if we have a grouped bar 
% chart, but this also works with only 2 or more groups, so we need to create 
% an extra data group
ganges_vuln_colors = ganges_vuln(1:3); 
ganges_vuln_colors(2,:) = zeros;
% in this case you'll have to manually try out the right locations of bars and 
% labels and the relative width of bars to make them look nice. Might be
% easier to do in Illustrator, if you need just the one chart!
figure 
barh([0.40,0.55], ganges_vuln_stacked, 0.70, 'stacked'); 
%0.4 and 0.55 refer to the vertical axis, they are the locations of the extra 
%stacked bar (0.4) and the actual stacked bar (0.55)
hold on %with hold on we are able to draw another plot in the same axes
barh(1:2, ganges_vuln_colors, 0.5, 'grouped'); 
%1 and 2 are the locations of the parameter values group and the extra
%group we created
% set the Ylim so that either of the extra data sets is not inside our plot
set(gca,'Ylim',[0.45 1.4],'Xlim',[0, 1],'YTick', [0.55, 0.7767, 1.0033, 1.23],'YTickLabel',new_label_vuln, 'Ticklength', [0,0]);


% and the same thing for resilience
ganges_res_stacked(1,1:3) = zeros;
ganges_res_stacked(2,:) = ganges_res(1:3) / 3;

ganges_res_colors = ganges_res(1:3); 
ganges_res_colors(2,:) = zeros;

new_label_res = ['total' title_res_parameters(1:3)];

figure 
barh([0.40,0.55], ganges_res_stacked, 0.70, 'stacked'); 

hold on 
barh(1:2, ganges_res_colors, 0.5, 'grouped'); 
set(gca,'Ylim',[0.45 1.4],'Xlim',[0, 1],'YTick', [0.55, 0.7767, 1.0033, 1.23],'YTickLabel',new_label_res, 'Ticklength', [0,0]);

%% EXAMPLE: Map the parameters and the vulnerability inside one basin

% open the georeference that matches with our parameters variable
load support_05dgr.mat R_05
load clip_zones.mat basin_names

[ganges_res_parameters_norm, ganges_R] = cut_basin_raster(res_parameters_norm, R_05, 'Ganges-Brahmaputra', basin_names, basins_05dgr);
[ganges_vuln_parameters_norm, ganges_R] = cut_basin_raster(vuln_parameters_norm, R_05, 'Ganges-Brahmaputra', basin_names, basins_05dgr);

ganges_res_index = nanmean(ganges_res_parameters_norm,3);
ganges_vuln_index = nanmean(ganges_vuln_parameters_norm,3);

% parameters
title_res = {'Governance','Economy (GDP/cap)','Social capital (HDI)'}; % NOTE: curly brackets needed for string type array
title_vuln = {'Human Footprint','Natural Hazards','Water Scarcity'}; % NOTE: curly brackets needed for string type array

figure
for i = 1:3
    subplot(2,3,i)
    imagesc(ganges_res_parameters_norm(:,:,i))
    title(title_res{i})
    colorbar
    subplot(2,3,i+3)
    imagesc(ganges_vuln_parameters_norm(:,:,i))
    title(title_vuln{i})
    colorbar
end


% index
figure;
    subplot(1,2,1)
    imagesc(ganges_res_index);
    colorbar
    title('resilience index');
    subplot(1,2,2)
    imagesc(ganges_vuln_index);
    colorbar
    title('vulnerability index');