clear
close all
clc
format compact
format short

set_list;

% Total PM 2.5
[~, ~, raw] = xlsread('plot.xlsx','PM_total','B2:F31');
data = reshape([raw{:}],size(raw));
pm_total_r_2010 = data(:,1);
pm_total_r_2030_bau = data(:,2);
pm_total_r_2030_p3 = data(:,3);
pm_total_r_2030_p4 = data(:,4);
pm_total_r_2030_p5 = data(:,5);
clearvars data raw;

pm_r_diff_p3_pctg = (pm_total_r_2030_p3 - pm_total_r_2030_bau)./pm_total_r_2030_bau;
pm_r_diff_p4_pctg = (pm_total_r_2030_p4 - pm_total_r_2030_bau)./pm_total_r_2030_bau;
pm_r_diff_p5_pctg = (pm_total_r_2030_p5 - pm_total_r_2030_bau)./pm_total_r_2030_bau;


%% Figure 3a: PM 2.5 in 2010
figure(1); clf; hold on; box on;
set(gcf, 'unit', 'inch', 'pos', [2.02    2.02    5.85    4.03]);
ax = worldmap('world');

latlim = [13 56];
lonlim = [70 140];
setm(ax, 'MapLatLimit', latlim);
setm(ax, 'MapLonLimit', lonlim);
setm(ax, 'MapProjection', 'miller');

load coast;
geoshow(flipud(lat), flipud(long), 'DisplayType', 'polygon', 'FaceColor', [0.8 0.97 1], 'edgecolor', [0.7 0.7 0.7]); %flipud lat and lon to fill the sea

load country_boundary;
plotm(country_boundary_lat, country_boundary_lon, 'color', [0.8 0.8 0.8]); % country line of whole world

clr = colormap(parula(100));
min_tick = 0;
max_tick = 120;

% Assign province color based on data
tmp = pm_total_r_2010;
tmp(tmp>max_tick) = max_tick;
tmp(tmp<min_tick) = min_tick;
scaling_color = [interp1(linspace(min_tick, max_tick, length(clr)), clr(:,1), tmp), ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,2), tmp), ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,3), tmp)];
scaling_color = [scaling_color;[0.9 0.9 0.9]]; % add 'XZ' as grey mannually

load china_province;
faceColors = makesymbolspec('Polygon',{'INDEX', [1 numel(china_province)], 'FaceColor', scaling_color});
geoshow(gca, china_province, 'DisplayType', 'polygon', 'SymbolSpec', faceColors, 'edgecolor', 'k')
framem off; gridm off; mlabel off; plabel off

set(gca, 'pos', [-0.16    -0.27    1.3    1.4]);
set(gcf, 'color', 'w');

% Lable the province names
latlim = getm(ax, 'MapLatLimit');
lonlim = getm(ax, 'MapLonLimit');
label_lat = [china_province.LabelLat];
label_lon = [china_province.LabelLon];
tf = ingeoquad(label_lat, label_lon, latlim, lonlim);
% textm(label_lat(tf(1:end-1)), label_lon(tf(1:end-1)), r1, 'HorizontalAlignment', 'center', 'fontsize', 8, 'fontweight', 'bold', 'color', 'k');
% textm(label_lat(end), label_lon(end), 'Tibet', 'HorizontalAlignment', 'center', 'fontsize', 8, 'fontweight', 'bold', 'color', [1 1 1]*0.8);

% Set colormap
color_tick = 50;
tmp = linspace(min_tick, max_tick, color_tick);
scaling_color = [interp1(linspace(min_tick, max_tick, length(clr)), clr(:,1), tmp)', ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,2), tmp)', ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,3), tmp)'];
colormap(scaling_color);  
caxis([min_tick, max_tick]);

% Set colorbar
c = colorbar;
c.Position = [0.82   0.1    0.03    0.35];
c.Limits = [min_tick, max_tick];
%c.Ticks = linspace(min_tick, max_tick, 6);
c.Ticks = [];
%c.Label.String = {'Pop.-Weighted PM_{2.5}', 'in 2010 (\mug/m_3)'};
%c.Label.FontSize = 13;

%textm(52, 77, '(a)', 'HorizontalAlignment', 'center', 'fontsize', 30, 'fontweight', 'bold', 'color', [0 0 0]);

export_fig fig3a -painters -png -r300;

%% Figure 2b: PM 2.5 in 2030 BAU
figure(2); clf; hold on; box on;
set(gcf, 'unit', 'inch', 'pos', [2.02    2.02    5.85    4.03]);
ax = worldmap('world');

latlim = [13 56];
lonlim = [70 140];
setm(ax, 'MapLatLimit', latlim);
setm(ax, 'MapLonLimit', lonlim);
setm(ax, 'MapProjection', 'miller');

load coast;
geoshow(flipud(lat), flipud(long), 'DisplayType', 'polygon', 'FaceColor', [0.8 0.97 1], 'edgecolor', [0.7 0.7 0.7]); %flipud lat and lon to fill the sea

load country_boundary;
plotm(country_boundary_lat, country_boundary_lon, 'color', [0.8 0.8 0.8]); % country line of whole world

% Assign province color based on data
tmp = pm_total_r_2030_bau;
tmp(tmp>max_tick) = max_tick;
tmp(tmp<min_tick) = min_tick;
scaling_color = [interp1(linspace(min_tick, max_tick, length(clr)), clr(:,1), tmp), ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,2), tmp), ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,3), tmp)];
scaling_color = [scaling_color;[0.9 0.9 0.9]]; % add 'XZ' as grey mannually

load china_province;
faceColors = makesymbolspec('Polygon',{'INDEX', [1 numel(china_province)], 'FaceColor', scaling_color});
geoshow(gca, china_province, 'DisplayType', 'polygon', 'SymbolSpec', faceColors, 'edgecolor', 'k')
framem off; gridm off; mlabel off; plabel off

set(gca, 'pos', [-0.16    -0.27    1.3    1.4]);
set(gcf, 'color', 'w');

% Lable the province names
latlim = getm(ax, 'MapLatLimit');
lonlim = getm(ax, 'MapLonLimit');
label_lat = [china_province.LabelLat];
label_lon = [china_province.LabelLon];
tf = ingeoquad(label_lat, label_lon, latlim, lonlim);
% textm(label_lat(tf(1:end-1)), label_lon(tf(1:end-1)), r1, 'HorizontalAlignment', 'center', 'fontsize', 8, 'fontweight', 'bold', 'color', 'k');
% textm(label_lat(end), label_lon(end), 'Tibet', 'HorizontalAlignment', 'center', 'fontsize', 8, 'fontweight', 'bold', 'color', [1 1 1]*0.8);

% Set colormap
color_tick = 50;
tmp = linspace(min_tick, max_tick, color_tick);
scaling_color = [interp1(linspace(min_tick, max_tick, length(clr)), clr(:,1), tmp)', ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,2), tmp)', ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,3), tmp)'];
colormap(scaling_color);  
caxis([min_tick, max_tick]);

% Set colorbar
c = colorbar;
c.Position = [0.82   0.1    0.03    0.35];
c.Limits = [min_tick, max_tick];
%c.Ticks = linspace(min_tick, max_tick, 6);
c.Ticks = [];
%c.Label.String = {'Pop.-Weighted PM_{2.5}', 'in 2030 No Policy (\mug/m_3)'};
%c.Label.FontSize = 13;

%textm(52, 77, '(b)', 'HorizontalAlignment', 'center', 'fontsize', 30, 'fontweight', 'bold', 'color', [0 0 0]);

export_fig fig3b -painters -png -r300;


%% Figure 2c: PM2.5 pctg difference in 3 percent policy
figure(3); clf; hold on; box on;
set(gcf, 'unit', 'inch', 'pos', [2.02    2.02    5.85    4.03]);
ax = worldmap('world');

latlim = [13 56];
lonlim = [70 140];
setm(ax, 'MapLatLimit', latlim);
setm(ax, 'MapLonLimit', lonlim);
setm(ax, 'MapProjection', 'miller');

load coast;
geoshow(flipud(lat), flipud(long), 'DisplayType', 'polygon', 'FaceColor', [0.8 0.97 1], 'edgecolor', [0.7 0.7 0.7]); %flipud lat and lon to fill the sea

load country_boundary;
plotm(country_boundary_lat, country_boundary_lon, 'color', [0.8 0.8 0.8]);

clr = [     0    0.2667    0.1059
            0    0.4275    0.1725
       0.1373    0.5451    0.2706
       0.2549    0.6706    0.3647
       0.4549    0.7686    0.4627
       0.6314    0.8510    0.6078
       0.7804    0.9137    0.7529
       0.8980    0.9608    0.8784
       0.9686    0.9882    0.9608
       1.0000    1.0000    1.0000];
min_tick = -0.25*100;
max_tick = 0;

% assign province color based on data
tmp = pm_r_diff_p3_pctg*100;
tmp(tmp>max_tick) = max_tick;
tmp(tmp<min_tick) = min_tick;
scaling_color = [interp1(linspace(min_tick, max_tick, length(clr)), clr(:,1), tmp), ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,2), tmp), ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,3), tmp)];
scaling_color = [scaling_color;[0.9 0.9 0.9]]; % add 'XZ' as grey mannually

load china_province;
faceColors = makesymbolspec('Polygon',{'INDEX', [1 numel(china_province)], 'FaceColor', scaling_color});
geoshow(gca, china_province, 'DisplayType', 'polygon', 'SymbolSpec', faceColors, 'edgecolor', 'k')
framem off; gridm off; mlabel off; plabel off

set(gca, 'pos', [-0.16    -0.27    1.3    1.4]);
set(gcf, 'color', 'w');

% Lable the province names
latlim = getm(ax, 'MapLatLimit');
lonlim = getm(ax, 'MapLonLimit');
label_lat = [china_province.LabelLat];
label_lon = [china_province.LabelLon];
tf = ingeoquad(label_lat, label_lon, latlim, lonlim);
% textm(label_lat(tf(1:end-1)), label_lon(tf(1:end-1)), r1, 'HorizontalAlignment', 'center', 'fontsize', 8, 'fontweight', 'bold', 'color', 'k');
% textm(label_lat(end), label_lon(end), 'Tibet', 'HorizontalAlignment', 'center', 'fontsize', 8, 'fontweight', 'bold', 'color', [1 1 1]*0.8);

% Set colormap
tmp = linspace(min_tick, max_tick, 50);
scaling_color = [interp1(linspace(min_tick, max_tick, length(clr)), clr(:,1), tmp)', ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,2), tmp)', ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,3), tmp)'];
colormap(scaling_color);  
caxis([min_tick, max_tick]);

% Set colorbar
c = colorbar;
c.Position = [0.82   0.1    0.03    0.35];
c.Limits = [min_tick, max_tick];
%c.Ticks = linspace(min_tick, max_tick, 6);
c.Ticks = [];
%c.Label.String = {'\DeltaPM_{2.5} in 2030 (%)', '(3% Policy vs. No Policy)'};
%c.Label.FontSize = 13;

%textm(52, 77, '(c)', 'HorizontalAlignment', 'center', 'fontsize', 30, 'fontweight', 'bold', 'color', [0 0 0]);
%textm(20, 75, 'National Change = -4.7%', 'HorizontalAlignment', 'left', 'fontsize', 18, 'fontweight', 'normal', 'color', [0 0 0]);

export_fig fig3c -painters -png -r300;


%% Figure 2d: PM2.5 pctg difference in 4 percent policy
figure(4); clf; hold on; box on;
set(gcf, 'unit', 'inch', 'pos', [2.02    2.02    5.85    4.03]);
ax = worldmap('world');

latlim = [13 56];
lonlim = [70 140];
setm(ax, 'MapLatLimit', latlim);
setm(ax, 'MapLonLimit', lonlim);
setm(ax, 'MapProjection', 'miller');

load coast;
geoshow(flipud(lat), flipud(long), 'DisplayType', 'polygon', 'FaceColor', [0.8 0.97 1], 'edgecolor', [0.7 0.7 0.7]); %flipud lat and lon to fill the sea

load country_boundary;
plotm(country_boundary_lat, country_boundary_lon, 'color', [0.8 0.8 0.8]);

% assign province color based on data
tmp = pm_r_diff_p4_pctg*100;
tmp(tmp>max_tick) = max_tick;
tmp(tmp<min_tick) = min_tick;
scaling_color = [interp1(linspace(min_tick, max_tick, length(clr)), clr(:,1), tmp), ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,2), tmp), ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,3), tmp)];
scaling_color = [scaling_color;[0.9 0.9 0.9]]; % add 'XZ' as grey mannually

load china_province;
faceColors = makesymbolspec('Polygon',{'INDEX', [1 numel(china_province)], 'FaceColor', scaling_color});
geoshow(gca, china_province, 'DisplayType', 'polygon', 'SymbolSpec', faceColors, 'edgecolor', 'k')
framem off; gridm off; mlabel off; plabel off

set(gca, 'pos', [-0.16    -0.27    1.3    1.4]);
set(gcf, 'color', 'w');

% Lable the province names
latlim = getm(ax, 'MapLatLimit');
lonlim = getm(ax, 'MapLonLimit');
label_lat = [china_province.LabelLat];
label_lon = [china_province.LabelLon];
tf = ingeoquad(label_lat, label_lon, latlim, lonlim);
% textm(label_lat(tf(1:end-1)), label_lon(tf(1:end-1)), r1, 'HorizontalAlignment', 'center', 'fontsize', 8, 'fontweight', 'bold', 'color', 'k');
% textm(label_lat(end), label_lon(end), 'Tibet', 'HorizontalAlignment', 'center', 'fontsize', 8, 'fontweight', 'bold', 'color', [1 1 1]*0.8);

% Set colormap
tmp = linspace(min_tick, max_tick, 50);
scaling_color = [interp1(linspace(min_tick, max_tick, length(clr)), clr(:,1), tmp)', ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,2), tmp)', ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,3), tmp)'];
colormap(scaling_color);  
caxis([min_tick, max_tick]);

% Set colorbar
c = colorbar;
c.Position = [0.82   0.1    0.03    0.35];
c.Limits = [min_tick, max_tick];
%c.Ticks = linspace(min_tick, max_tick, 6);
c.Ticks = [];
%c.Label.String = {'\DeltaPM_{2.5} in 2030 (%)', '(4% Policy vs. No Policy)'};
%c.Label.FontSize = 13;

%textm(52, 77, '(d)', 'HorizontalAlignment', 'center', 'fontsize', 30, 'fontweight', 'bold', 'color', [0 0 0]);
%textm(20, 75, 'National Change = -11.9%', 'HorizontalAlignment', 'left', 'fontsize', 18, 'fontweight', 'normal', 'color', [0 0 0]);

export_fig fig3d -painters -png -r300;


%% Figure 2e: PM2.5 pctg difference in 4 percent policy
figure(5); clf; hold on; box on;
set(gcf, 'unit', 'inch', 'pos', [2.02    2.02    5.85    4.03]);
ax = worldmap('world');

latlim = [13 56];
lonlim = [70 140];
setm(ax, 'MapLatLimit', latlim);
setm(ax, 'MapLonLimit', lonlim);
setm(ax, 'MapProjection', 'miller');

load coast;
geoshow(flipud(lat), flipud(long), 'DisplayType', 'polygon', 'FaceColor', [0.8 0.97 1], 'edgecolor', [0.7 0.7 0.7]); %flipud lat and lon to fill the sea

load country_boundary;
plotm(country_boundary_lat, country_boundary_lon, 'color', [0.8 0.8 0.8]);

% assign province color based on data
tmp = pm_r_diff_p5_pctg*100;
tmp(tmp>max_tick) = max_tick;
tmp(tmp<min_tick) = min_tick;
scaling_color = [interp1(linspace(min_tick, max_tick, length(clr)), clr(:,1), tmp), ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,2), tmp), ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,3), tmp)];
scaling_color = [scaling_color;[0.9 0.9 0.9]]; % add 'XZ' as grey mannually

load china_province;
faceColors = makesymbolspec('Polygon',{'INDEX', [1 numel(china_province)], 'FaceColor', scaling_color});
geoshow(gca, china_province, 'DisplayType', 'polygon', 'SymbolSpec', faceColors, 'edgecolor', 'k')
framem off; gridm off; mlabel off; plabel off

set(gca, 'pos', [-0.16    -0.27    1.3    1.4]);
set(gcf, 'color', 'w');

% Lable the province names
latlim = getm(ax, 'MapLatLimit');
lonlim = getm(ax, 'MapLonLimit');
label_lat = [china_province.LabelLat];
label_lon = [china_province.LabelLon];
tf = ingeoquad(label_lat, label_lon, latlim, lonlim);
% textm(label_lat(tf(1:end-1)), label_lon(tf(1:end-1)), r1, 'HorizontalAlignment', 'center', 'fontsize', 8, 'fontweight', 'bold', 'color', 'k');
% textm(label_lat(end), label_lon(end), 'Tibet', 'HorizontalAlignment', 'center', 'fontsize', 8, 'fontweight', 'bold', 'color', [1 1 1]*0.8);

% Set colormap
tmp = linspace(min_tick, max_tick, 50);
scaling_color = [interp1(linspace(min_tick, max_tick, length(clr)), clr(:,1), tmp)', ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,2), tmp)', ...
                 interp1(linspace(min_tick, max_tick, length(clr)), clr(:,3), tmp)'];
colormap(scaling_color);  
caxis([min_tick, max_tick]);

% Set colorbar
c = colorbar;
c.Position = [0.82   0.1    0.03    0.35];
c.Limits = [min_tick, max_tick];
%c.Ticks = linspace(min_tick, max_tick, 6);
c.Ticks = [];
%c.Label.String = {'\DeltaPM_{2.5} in 2030 (%)', '(5% Policy vs. No Policy)'};
%c.Label.FontSize = 13;

%textm(52, 77, '(e)', 'HorizontalAlignment', 'center', 'fontsize', 30, 'fontweight', 'bold', 'color', [0 0 0]);
%textm(20, 75, 'National Change = -18.7%', 'HorizontalAlignment', 'left', 'fontsize', 18, 'fontweight', 'normal', 'color', [0 0 0]);
export_fig fig3e -painters -png -r300;
