%% Lecture 6 - Socio-economic indicators


clc
clear
close all
addpath(genpath('path to your course folder here')); 

cd('path to your working directory here');

%
%% Load data
% water_scarcity_05dgr - Combined Water Scarcity Index for year 2015  (water shortage, WCI & water stress, WSI combined) based on Kummu et al (2016) and Wada et al (2014)
% gdp_per_capita_05dgr - GDP (US$) per capita PPP for year 2015  (Kummu et al 2018 - data available at https://doi.org/10.5061/dryad.dk1j0)
% hum_ftprint_05dgr - Human Footprint Index, year 2015 the relative human influence in terrestrial biomes  (0 = least influenced parts, 100 = most influenced parts of the biome) (Venter et al 2016, data available at http://doi.org/10.1038/sdata.2016.67)
% natural_haz_05dgr - global multihazard frequency and distribution average over longer time period (http://sedac.ciesin.columbia.edu/data/set/ndh-multihazard-frequency-distribution)
% hdi_05dgr - Human Development Index for year 2015 (Kummu et al 2018 - data available at https://doi.org/10.5061/dryad.dk1j0)
% Load data

load raster_input_week6.mat

%
%% Investigate what the data looks like
% 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 Scarcity (combined water stress and shortage)')

subplot(2,3,5);
imagesc(hdi_05dgr);
set(gca, 'CLim', [0, 2]);
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); 

% For region 13 (ice), global average

temp(13,1) = nanmean(governance_indicator(:,3));

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 j

%
%% Governance data to raster
% Load country raster

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

% 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
% There are some negative values in the GDP dataset, (-9 for areas without data) and thus, we need get rid of those first

gdp_per_capita_05dgr(gdp_per_capita_05dgr== -9) = NaN;
gdp_per_capita_05dgr_log = log(gdp_per_capita_05dgr); 

%
%% Vulnerability and adaptive capacity index for each grid cell
% Adaptive capacity index is calculated based on 3 societal indicators:
% Governance
% GDP/cap
% HDI
% Create an empty 3-layered variable

cap_parameters = zeros(size(adm0_05dgr,1),size(adm0_05dgr,2),3);

% A. collect the indicators:

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

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

% 3. social capital: HDI
cap_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

cap_parameters_norm = zeros(size(cap_parameters));
vuln_parameters_norm = zeros(size(vuln_parameters));

for i = 1:3
    temp_cap = cap_parameters(:,:,i);
    temp_vuln = vuln_parameters(:,:,i);
   
    % normalise adaptive capacity index in question
    
    temp_5perc = prctile(temp_cap(:),5);
    temp_95perc = prctile(temp_cap(:),95);
    cap_index = (temp_cap - temp_5perc) ./ (temp_95perc-temp_5perc);
      
    % normalise vulnerability index
    % note: in case water scarcity, no normalisation needed as scale is between 0 and 1
    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
    
    cap_index(cap_index > 1) = 1;
    cap_index(cap_index < 0) = 0;
    vuln_index(vuln_index > 1) = 1;
    vuln_index(vuln_index < 0) = 0;
    
    
    cap_parameters_norm(:,:,i) = cap_index;
    vuln_parameters_norm(:,:,i) = vuln_index;

end

% Calculate vulnerability and adaptive capacity indices

vuln_index = nanmean(vuln_parameters_norm,3);
cap_index = nanmean(cap_parameters_norm,3);

%
%% Plot the parameters and the vulnerability index
% Parameters

title_cap_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(cap_parameters_norm(:,:,i))
    title(title_cap_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(cap_index);
    colorbar
    title('adaptive capacity index');
    subplot(1,2,2)
    imagesc(vuln_index);
    colorbar
    title('vulnerability index');

%
%% Calculate the individual vulnerability and adaptive capacity parameters for one basin
% For one basin, e.g. Nile (id = 19)
% Load basins

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

% Adaptive capacity 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');

% Select basin using a logical variable

basin_logical = basins_05dgr == 19;

% For adaptive capacity parameters calculate the sum of population weighted parameters inside nile basin and divide with total population of the basin

for i = 1:3
    cap_parameters_norm_weighted(:,:,i) = cap_parameters_norm(:,:,i) .* total_pop_sel(:,:,5);
    temp_weighted = basin_logical .* cap_parameters_norm_weighted(:,:,i);  
    temp_pop = basin_logical .* total_pop_sel(:,:,5);
    nile_cap(1,i) = nansum(nansum(temp_weighted))./nansum(nansum(temp_pop));
end

% For vulnerability parameters calculate the sum of area weighted parameters inside nile 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;
    nile_vuln(1,i) = nansum(nansum(temp_weighted))./nansum(nansum(temp_area));
end
clear temp*

nile_cap(1,4) = nanmean(nile_cap);    
nile_vuln(1,4) = nanmean(nile_vuln);

%
%% EXAMPLE: Plot the results as a bar chart
% Bar chart of different parameters

figure
title_vuln = [title_vuln_parameters,'total'];
barh(nile_vuln, 0.5); % 0.5 is the relative bar width 
set(gca,'YTickLabel',title_vuln);

% Stacked bar chart of vuln index
% Stacked bar chart only works with 2 or more stacked bars, so we need to fool matlab by creating an extra dataset:

nile_vuln_stacked(1,1:3) = zeros; 
nile_vuln_stacked(2,:) = nile_vuln(1:3) / 3; 

figure
h = barh(nile_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], nile_vuln_stacked, 0.5, 'stacked');
hold on
h = barh(3:5, nile_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

nile_vuln_colors = nile_vuln(1:3); 
nile_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!
% 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)

figure 
h1 = barh([0.40,0.55], nile_vuln_stacked, 0.70, 'stacked'); 

% Here, 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

hold on %with hold on we are able to draw another plot in the same axes
h2 = barh(1:2, nile_vuln_colors, 0.5, 'grouped'); 
h2(1).FaceColor = h1(1).FaceColor;
h2(2).FaceColor = h1(2).FaceColor;
h2(3).FaceColor = h1(3).FaceColor;
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]);

% The same thing for adaptive capacity

nile_cap_stacked(1,1:3) = zeros;
nile_cap_stacked(2,:) = nile_cap(1:3) / 3;

nile_cap_colors = nile_cap(1:3); 
nile_cap_colors(2,:) = zeros;

new_label_cap = ['total' title_cap_parameters(1:3)];

figure 
h3 = barh([0.40,0.55], nile_cap_stacked, 0.70, 'stacked'); 

hold on 
h4 = barh(1:2, nile_cap_colors, 0.5, 'grouped'); 
h4(1).FaceColor = h3(1).FaceColor;
h4(2).FaceColor = h3(2).FaceColor;
h4(3).FaceColor = h3(3).FaceColor;
set(gca,'Ylim',[0.45 1.4],'Xlim',[0, 1],'YTick', [0.55, 0.7767, 1.0033, 1.23],'YTickLabel',new_label_cap, '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

% Cut the adaptive capacity and vulnerability to one basin

[nile_cap_parameters_norm, nile_R] = cut_basin_raster(cap_parameters_norm, R_05, 'Nile', basin_names, basins_05dgr);
[nile_vuln_parameters_norm, nile_R] = cut_basin_raster(vuln_parameters_norm, R_05, 'Nile', basin_names, basins_05dgr);

nile_cap_index = nanmean(nile_cap_parameters_norm,3);
nile_vuln_index = nanmean(nile_vuln_parameters_norm,3);

% Parameters

title_cap = {'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(nile_cap_parameters_norm(:,:,i))
    title(title_cap{i})
    colorbar
    subplot(2,3,i+3)
    imagesc(nile_vuln_parameters_norm(:,:,i))
    title(title_vuln{i})
    colorbar
end

% index

figure;
    subplot(1,2,1)
    imagesc(nile_cap_index);
    colorbar
    title('adaptive capacity index');
    subplot(1,2,2)
    imagesc(nile_vuln_index);
    colorbar
    title('vulnerability index');