%Read in the sample data
nc=ncgeodataset('GDP_PPP_30arcsec_v3.nc');
dirvar=nc.geovariable('pr')
 %dir=dirvar.data(:,:,:);
g=dirvar.grid_interop(:,:);
timebnd=nc.geovariable('time_bnds');%specific for UKESM1 use time_bnds for other cases.

%Manage the location
UKESM1_Cord=nan(length(g.lon)*length(g.lat), 4);
k=0;
for i=1:length(g.lon)
    for j=1:length(g.lat)
        k=k+1;
        UKESM1_Cord(k,1)=g.lon(i)-360;
        UKESM1_Cord(k,2)=g.lat(j);
        UKESM1_Cord(k,3)=i;
        UKESM1_Cord(k,4)=j;
    end
end


nc=ncgeodataset('GDP_PPP_30arcsec_v3.nc');
var=nc.geovariable('GDP_PPP');
g=var.grid_interop(1,:,:);
glonrange=find(g.lon>=min(min(trmm_pos(:,:,1))) & g.lon<=max(max(trmm_pos(:,:,1))));
glatrange=find(g.lat>=min(min(trmm_pos(:,:,2))) & g.lat<=max(max(trmm_pos(:,:,2))));
USGDP90=nan(4140,6900);
USGDP90(:,:)=var.data(1,glatrange,glonrange);
USGDP15=nan(4140,6900);
USGDP15(:,:)=var.data(3,glatrange,glonrange);
glon(:,1)=g.lon(glonrange);
glat(:,1)=g.lat(glatrange);

GDPGRID=nan(4140,6900,2);
for i=1:4140
    for j=1:6900
    GDPGRID(i,j,1)=glon(j,1);
    GDPGRID(i,j,2)=glat(i,1);
    end
end
 

GDP1990=NaN(139,231);
GDP2015=NaN(139,231);
for i=1:139
    for j=1:231
        if Mask_G(i,j)==1
[ggdp(:,1),ggdp(:,2)]=find(GDPGRID(:,:,1)>=(trmm_pos(i,j,1)-0.25)&GDPGRID(:,:,1)<=trmm_pos(i,j,1)...
    &GDPGRID(:,:,2)>=trmm_pos(i,j,2)&GDPGRID(:,:,2)<=(trmm_pos(i,j,2)+0.25));
GDP1990(i,j)=nansum(nansum(diag(USGDP90(ggdp(:,1),ggdp(:,2)))));
GDP2015(i,j)=nansum(nansum(diag(USGDP15(ggdp(:,1),ggdp(:,2)))));
clear ggdp
        end
    end
    i
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 subplot(1,2,1)
m_proj('Lambert','lon',NEWLONLIMS,'lat',NEWLATLIMS)
m_pcolor(trmm_pos(:,:,1), trmm_pos(:,:,2),log(GDP1990)./EventRisk100_NP1);shading flat;
m_grid('linewi',2,'tickdir','out')
% hold on
% m_plot(bndlon4(:)-0.125,bndlat4(:)+0.125,'x','color','k');
caxis([0 8.1943])
title('a. Early TCP Risk (> 100mm) Exposure to GDP')
for k=1:length(M.ncst)
m_line(M.ncst{k}(:,1),M.ncst{k}(:,2),'color','k','linewidth', 0.5);
end
colormap(othercolor('RdPu9'))
set(gca,'fontsize',14);


subplot(1,2,2)
m_proj('Lambert','lon',NEWLONLIMS,'lat',NEWLATLIMS)
m_pcolor(trmm_pos(:,:,1), trmm_pos(:,:,2),log(GDP2015)./EventRisk100_NP2);shading flat;
m_grid('linewi',2,'tickdir','out')
% hold on
% m_plot(bndlon4(:)-0.125,bndlat4(:)+0.125,'x','color','k');
caxis([0 8.1943])
title('b. Late TCP Risk (> 100mm) Exposure to GDP')
for k=1:length(M.ncst)
m_line(M.ncst{k}(:,1),M.ncst{k}(:,2),'color','k','linewidth', 0.5);
end
colormap(othercolor('RdPu9'))
set(gca,'fontsize',14);


set(gca,'fontsize',14);
h = colorbar;
colorTitleHandle = get(h,'Title');
titleString = 'TCP Risk Exposure';
set(colorTitleHandle ,'String',titleString);
 set(h, 'Position', [.915 .11 .0181 .8550]);
 set(gca,'fontsize',14);

 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
 
 POPEI2020=log(P2020)./EventRisk100_NP2;
 POPEI1970=log(P1970)./EventRisk100_NP1;
 
 GDPEI2015=log(GDP2015)./EventRisk100_NP2;
 GDPEI1990=log(GDP1990)./EventRisk100_NP1;
 
 POPEI22020=log(P2020)./EventRisk200_NP2;
 POPEI21970=log(P1970)./EventRisk200_NP1;
 
 GDPEI22015=log(GDP2015)./EventRisk200_NP2;
 GDPEI21990=log(GDP1990)./EventRisk200_NP1;
 
 
 subplot(1,2,1)
 PopTCPEI=(log(P2020)./EventRisk100_NP2-log(P1970)./EventRisk100_NP1);
 PopTCPEI(PopTCPEI==Inf)=NaN;
m_proj('Lambert','lon',NEWLONLIMS,'lat',NEWLATLIMS)
m_pcolor(trmm_pos(:,:,1), trmm_pos(:,:,2),PopTCPEI);shading flat;
m_grid('linewi',2,'tickdir','out')
% hold on
% m_plot(bndlon4(:)-0.125,bndlat4(:)+0.125,'x','color','k');
caxis([-2.2 2.2])
title('a. Change in Population Exposure to TCP Risk (> 100mm)')
for k=1:length(M.ncst)
m_line(M.ncst{k}(:,1),M.ncst{k}(:,2),'color','k','linewidth', 0.5);
end
colormap(jet)
set(gca,'fontsize',14);
h=colorbar;
colorTitleHandle = get(h,'Title');
titleString = 'TCP Risk Exposure';

subplot(1,2,2)
GDPTCPEI=(log(GDP2015)./EventRisk100_NP2-log(GDP1990)./EventRisk100_NP1);
GDPTCPEI(GDPTCPEI==Inf)=NaN;

m_proj('Lambert','lon',NEWLONLIMS,'lat',NEWLATLIMS)
m_pcolor(trmm_pos(:,:,1), trmm_pos(:,:,2),GDPTCPEI);shading flat;
m_grid('linewi',2,'tickdir','out')
% hold on
% m_plot(bndlon4(:)-0.125,bndlat4(:)+0.125,'x','color','k');
caxis([-4 4])
title('b. Change in GDP Exposure to TCP Risk (> 100mm)')
for k=1:length(M.ncst)
m_line(M.ncst{k}(:,1),M.ncst{k}(:,2),'color','k','linewidth', 0.5);
end
colormap(jet)
set(gca,'fontsize',14);
h=colorbar;
colorTitleHandle = get(h,'Title');
titleString = 'TCP Risk Exposure';

% set(gca,'fontsize',14);
% h = colorbar;
% colorTitleHandle = get(h,'Title');
% titleString = 'TCP Risk Exposure';
% set(colorTitleHandle ,'String',titleString);
%  set(h, 'Position', [.915 .11 .0181 .8550]);
%  set(gca,'fontsize',14);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 subplot(1,2,1)
 PopTCPEIP=(log(P2020)./EventRisk100_NP2-log(P1970)./EventRisk100_NP1)./(log(P1970)./EventRisk100_NP1).*100;
 PopTCPEIP(PopTCPEIP==Inf)=NaN;
m_proj('Lambert','lon',NEWLONLIMS,'lat',NEWLATLIMS)
m_pcolor(trmm_pos(:,:,1), trmm_pos(:,:,2),PopTCPEIP);shading flat;
m_grid('linewi',2,'tickdir','out')
% hold on
% m_plot(bndlon4(:)-0.125,bndlat4(:)+0.125,'x','color','k');
%caxis([-2.2 2.2])
title('a. Change in Population Exposure to TCP Risk (> 100mm)')
for k=1:length(M.ncst)
m_line(M.ncst{k}(:,1),M.ncst{k}(:,2),'color','k','linewidth', 0.5);
end
colormap(jet)
set(gca,'fontsize',14);
h=colorbar;
colorTitleHandle = get(h,'Title');
titleString = 'TCP Risk Exposure';

subplot(1,2,2)
GDPTCPEIP=(log(GDP2015)./EventRisk100_NP2-log(GDP1990)./EventRisk100_NP1)/(log(GDP1990)./EventRisk100_NP1);
GDPTCPEIP(GDPTCPEIP==Inf)=NaN;

m_proj('Lambert','lon',NEWLONLIMS,'lat',NEWLATLIMS)
m_pcolor(trmm_pos(:,:,1), trmm_pos(:,:,2),GDPTCPEIP);shading flat;
m_grid('linewi',2,'tickdir','out')
% hold on
% m_plot(bndlon4(:)-0.125,bndlat4(:)+0.125,'x','color','k');
% caxis([-4 4])
title('b. Change in GDP Exposure to TCP Risk (> 100mm)')
for k=1:length(M.ncst)
m_line(M.ncst{k}(:,1),M.ncst{k}(:,2),'color','k','linewidth', 0.5);
end
colormap(jet)
set(gca,'fontsize',14);
h=colorbar;
colorTitleHandle = get(h,'Title');
titleString = 'TCP Risk Exposure';




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Pop_diff=log(P2020)-log(P1970);
subplot(1,2,1)
m_proj('Lambert','lon',NEWLONLIMS,'lat',NEWLATLIMS)
m_pcolor(trmm_pos(:,:,1), trmm_pos(:,:,2),Pop_diff);shading flat;
m_grid('linewi',2,'tickdir','out')
% hold on
% m_plot(bndlon4(:)-0.125,bndlat4(:)+0.125,'x','color','k');
caxis([-2.4 2.4])
for k=1:length(M.ncst)
m_line(M.ncst{k}(:,1),M.ncst{k}(:,2),'color','k','linewidth', 0.5);
end
title('a. Change in Population')
colormap(gca,'jet')
set(gca,'fontsize',14);
h1=colorbar;
colorTitleHandle = get(h1,'Title');
titleString = 'log(Population)';
set(colorTitleHandle ,'String',titleString);


GDP_diff=log(GDP2015)-log(GDP1990.*1.81);
subplot(1,2,2)
m_proj('Lambert','lon',NEWLONLIMS,'lat',NEWLATLIMS)
m_pcolor(trmm_pos(:,:,1), trmm_pos(:,:,2),GDP_diff);shading flat;
m_grid('linewi',2,'tickdir','out')
% hold on
% m_plot(bndlon4(:)-0.125,bndlat4(:)+0.125,'x','color','k');
caxis([0 1.5])
for k=1:length(M.ncst)
m_line(M.ncst{k}(:,1),M.ncst{k}(:,2),'color','k','linewidth', 0.5);
end
title('b. Change in GDP')
colormap(gca,flip(hot))
set(gca,'fontsize',14);
h2=colorbar;
colorTitleHandle = get(h2,'Title');
titleString = 'log(GDP)';
set(colorTitleHandle ,'String',titleString);


% m_proj('Lambert','lon',NEWLONLIMS,'lat',NEWLATLIMS)
% m_pcolor(trmm_pos(:,:,1), trmm_pos(:,:,2),log(GDP2015-GDP1990));shading flat;
% m_grid('linewi',2,'tickdir','out')
% % hold on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%More code for the comparison of state
for i=1:length(states)
    state_id{i} = find(trmmstate==states(i));
end

GDPTCPEIP_R=reshape(GDPTCPEIP,[32109,1]);
PopTCPEIP_R=reshape(PopTCPEIP,[32109,1]);

GDPTCPEI_R=reshape(GDPTCPEI,[32109,1]);
PopTCPEI_R=reshape(PopTCPEI,[32109,1]);

POPEI2020_R=reshape(POPEI2020,[32109,1]);
POPEI1970_R=reshape(POPEI1970,[32109,1]);
 
GDPEI2015_R=reshape(GDPEI2015,[32109,1]);
GDPEI1990_R=reshape(GDPEI1990,[32109,1]);

POPEI22020_R=reshape(POPEI22020,[32109,1]);
POPEI21970_R=reshape(POPEI21970,[32109,1]);
 
GDPEI22015_R=reshape(GDPEI22015,[32109,1]);
GDPEI21990_R=reshape(GDPEI21990,[32109,1]);



for i=1:49
POPTCP_EI(i,1)=nanmedian(PopTCPEI_R(state_id{i}));
POPTCP_EI(i,2)=nanmean(PopTCPEI_R(state_id{i}));
POPTCP_EI(i,3)=nanmax(PopTCPEI_R(state_id{i}));
POPTCP_EI(i,4)=nanmin(PopTCPEI_R(state_id{i}));
POPTCP_EI(i,5)=nanstd(PopTCPEI_R(state_id{i}));
POPTCP_EI(i,6)=sum(~isnan(PopTCPEI_R(state_id{i})));
GDPTCP_EI(i,1)=nanmedian(GDPTCPEI_R(state_id{i}));
GDPTCP_EI(i,2)=nanmean(GDPTCPEI_R(state_id{i}));
GDPTCP_EI(i,3)=nanmax(GDPTCPEI_R(state_id{i}));
GDPTCP_EI(i,4)=nanmin(GDPTCPEI_R(state_id{i}));
GDPTCP_EI(i,5)=nanstd(GDPTCPEI_R(state_id{i}));
GDPTCP_EI(i,6)=sum(~isnan(GDPTCPEI_R(state_id{i})));
end

for i=1:49
POPTCP_PC(i,1)=nanmedian(PopTCPEIP_R(state_id{i}));
POPTCP_PC(i,2)=nanmean(PopTCPEIP_R(state_id{i}));
POPTCP_PC(i,3)=nanmax(PopTCPEIP_R(state_id{i}));
POPTCP_PC(i,4)=nanmin(PopTCPEIP_R(state_id{i}));
POPTCP_PC(i,5)=nanstd(PopTCPEIP_R(state_id{i}));
POPTCP_PC(i,6)=sum(~isnan(PopTCPEIP_R(state_id{i})));
GDPTCP_PC(i,1)=nanmedian(GDPTCPEIP_R(state_id{i}));
GDPTCP_PC(i,2)=nanmean(GDPTCPEIP_R(state_id{i}));
GDPTCP_PC(i,3)=nanmax(GDPTCPEIP_R(state_id{i}));
GDPTCP_PC(i,4)=nanmin(GDPTCPEIP_R(state_id{i}));
GDPTCP_PC(i,5)=nanstd(GDPTCPEIP_R(state_id{i}));
GDPTCP_PC(i,6)=sum(~isnan(GDPTCPEIP_R(state_id{i})));
end

StateSerial=find(~isnan(POPTCP_PC(:,1)));
StateSerial2=find(~isnan(GDPTCP_PC2(:,1)));

POPTCP_Sel=POPTCP_PC(StateSerial,:);
GDPTCP_Sel=GDPTCP_PC(StateSerial,:);


POPTCPE_Sel=POPTCP_EI(StateSerial,:);
GDPTCPE_Sel=GDPTCP_EI(StateSerial,:);

EXISTPOP=nan(1038,21);
EXTSTGDP=nan(1038,21);




for i=1:49
POPTCP_PC(i,1)=nanmedian(PopTCPEIP_R(state_id{i}));
POPTCP_PC(i,2)=nanmean(PopTCPEIP_R(state_id{i}));
POPTCP_PC(i,3)=nanmax(PopTCPEIP_R(state_id{i}));
POPTCP_PC(i,4)=nanmin(PopTCPEIP_R(state_id{i}));
POPTCP_PC(i,5)=nanstd(PopTCPEIP_R(state_id{i}));
POPTCP_PC(i,6)=sum(~isnan(PopTCPEIP_R(state_id{i})));
GDPTCP_PC(i,1)=nanmedian(GDPTCPEIP_R(state_id{i}));
GDPTCP_PC(i,2)=nanmean(GDPTCPEIP_R(state_id{i}));
GDPTCP_PC(i,3)=nanmax(GDPTCPEIP_R(state_id{i}));
GDPTCP_PC(i,4)=nanmin(GDPTCPEIP_R(state_id{i}));
GDPTCP_PC(i,5)=nanstd(GDPTCPEIP_R(state_id{i}));
GDPTCP_PC(i,6)=sum(~isnan(GDPTCPEIP_R(state_id{i})));
end




%%%%%%%clear;clc;
% Create example data
A=rand(100,10);
B=rand(200,10);
C=rand(150,10);

% prepare data
POPT=cell(20,2);
for ii=1:20
Ac{ii}=EXISTPOP1970(:,ii+1);
Bc{ii}=EXISTPOP2020(:,ii+1);
end
POPT=vertcat(Ac,Bc);

xlab=states(StateSerial(2:21));
col=[102,255,255,200;
51,153,255, 200;];
col=col/255

multiple_boxplot(POPT',xlab,{'Early', 'Late'},col')
%title('Here it is!')

GDPT=cell(20,2);
for ii=1:20
Cc{ii}=EXISTGDP1990(:,ii+1);
Dc{ii}=EXISTGDP2015(:,ii+1);
end
GDPT=vertcat(Cc,Dc);

% xlab=states(StateSerial(2:21));
% col=[102,255,255,200;
% 51,153,255, 200;];
% col=col/255;

multiple_boxplot(GDPT',xlab,{'Early', 'Late'},col')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,1,1)
multiple_boxplot(POPT',xlab,{'Early', 'Late'},col');
ylabel('Exposure Index');
title('a.Change in TCP Exposure Index based on Population');
ylim([0 3.5])
set(gca,'fontsize',12);

subplot(2,1,2)
multiple_boxplot(GDPT',xlab,{'Early', 'Late'},col')
ylabel('Exposure Index');
title('b.TCP Exposure Index based on Gross Domestic Product (GDP)');
ylim([0 6.5])
set(gca,'fontsize',12);





bar(POPTCP_PC(~isnan(POPTCP_PC)))

subplot(1,2,1)
boxplot(EXISTPOP(:,2:21),states(StateSerial(2:21)))
ylim([-100 700])
yline(0,'--')

subplot(1,2,2)
boxplot(EXTSTGDP(:,2:21),states(StateSerial2(2:21)))
ylim([-100 700])
yline(0,'--')


% prepare data
POPGDPT=cell(20,2);
for ii=1:20
Ec{ii}=EXISTPOP(:,ii+1);
Fc{ii}=EXTSTGDP(:,ii+1);
end
POPGDPT=vertcat(Ec,Fc);

xlab=states(StateSerial(2:21));
col=[102,255,255, 200;
51,153,255, 200]
 col=col/255


% legend('no change')
 
 subplot(2,2,1)
 PopTCPEI=(log(P2020)./EventRisk100_NP2-log(P1970)./EventRisk100_NP1);
 PopTCPEI(PopTCPEI==Inf)=NaN;
m_proj('Lambert','lon',NEWLONLIMS,'lat',NEWLATLIMS)
m_pcolor(trmm_pos(:,:,1), trmm_pos(:,:,2),PopTCPEI);shading flat;
m_grid('linewi',2,'tickdir','out')
% hold on
% m_plot(bndlon4(:)-0.125,bndlat4(:)+0.125,'x','color','k');
caxis([-2.2 2.2])
title('a. Change in Population Exposure to TCP Risk (> 100mm)')
for k=1:length(M.ncst)
m_line(M.ncst{k}(:,1),M.ncst{k}(:,2),'color','k','linewidth', 0.5);
end
colormap(jet)
set(gca,'fontsize',10);
h=colorbar;
colorTitleHandle = get(h,'Title');
titleString = 'TCP Risk Exposure';

subplot(2,2,2)
GDPTCPEI=(log(GDP2015)./EventRisk100_NP2-log(GDP1990)./EventRisk100_NP1);
GDPTCPEI(GDPTCPEI==Inf)=NaN;

m_proj('Lambert','lon',NEWLONLIMS,'lat',NEWLATLIMS)
m_pcolor(trmm_pos(:,:,1), trmm_pos(:,:,2),GDPTCPEI);shading flat;
m_grid('linewi',2,'tickdir','out')
% hold on
% m_plot(bndlon4(:)-0.125,bndlat4(:)+0.125,'x','color','k');
caxis([-4 4])
title('b. Change in GDP Exposure to TCP Risk (> 100mm)')
for k=1:length(M.ncst)
m_line(M.ncst{k}(:,1),M.ncst{k}(:,2),'color','k','linewidth', 0.5);
end
colormap(jet)
set(gca,'fontsize',14);
h=colorbar;
colorTitleHandle = get(h,'Title');
titleString = 'TCP Risk Exposure';

subplot(2,2,3:4)
multiple_boxplot(POPGDPT',xlab,{'Population', 'GDP'},col')
ylim([-100 700])
yline(0,'--','HandleVisibility','off')
ylabel('Relative Change from early to late period (%)')