
addpath('../')


%===========================================================================
% LOAD CLIMATE INDICES

% Load full ensembles for the nino and amo variables
efl='cesm_lme010_r12000_p998_state_2ozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:24:35.mat';
% LOAD JJA FOR AMO ANALYSIS
%efl='cesm_lme010_r12000_p998_state_2ozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:11:15.mat';

pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';
%load([pth,efl],'r_o','r_f','state_tp','xmeta','Xa_m','Xa_ens','reconYrs','p_yrs')
load([pth,efl])


% What years to look at?
analyrs1=1000:1800;
analyrs2=1801:1925; % reference
analyrs3=1000:2000;
%analyrs3=800:2000;
%analyrs4=1000:1925;
analyrs4=1000:1800;
[~,~,ir1] = intersect(analyrs1,r_o:r_f); % indices of recon segment
[~,~,ir2] = intersect(analyrs2,r_o:r_f); % indices of recon segment
[~,~,ir3] = intersect(analyrs3,r_o:r_f); % indices of recon segment
[~,~,ir4] = intersect(analyrs4,r_o:r_f); % indices of recon segment
% monthly indices
ovrlprng=rng_intr([r_o,r_f],[analyrs1(1),analyrs1(end)]);
mon_r=r_o:(1/12):(r_f+11/12);
irm1_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
irm1_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;

ovrlprng=rng_intr([r_o,r_f],[analyrs2(1),analyrs2(end)]);
mon_r=r_o:(1/12):(r_f+11/12);
irm2_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
irm2_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;

ovrlprng=rng_intr([r_o,r_f],[analyrs3(1),analyrs3(end)]);
mon_r=r_o:(1/12):(r_f+11/12);
irm3_a=find(abs(ovrlprng(1)-mon_r)<1e-6);
irm3_b=find(abs(ovrlprng(2)-1/12-mon_r)<1e-6)+12;



k=find(state_tp=='o');
xa_amo_ens=permute(squeeze(Xa_ens(xmeta{k}.id_X(1),:,:)),[2 1]);
% Ensemble mean of AMO
xa_amo_mn=mean(xa_amo_ens,2);
amo_off=mean(mean(xa_amo_ens(ir2,:),2));


k=find(state_tp=='2');
xa_gmt_ens=permute(squeeze(Xa_ens(xmeta{k}.id_X(1),:,:)),[2 1]);

k=find(state_tp=='n'); % specific nino variable
if strcmp(xmeta{k}.indinfo{3}.indxnm,'Nino 3.4'); i=3; else; error('Wrong Nino specification'); end
xa_nino_ens=reshape(permute(Xa_ens((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:,:),[1 3 2]),12*reconYrs,length(p_yrs));
xa_nino_mn=reshape(Xa_m((xmeta{k}.id_X(1)+(i-1)*12):(xmeta{k}.id_X(1)+i*12-1),:),12*reconYrs,1);
% Nino offset with respect to reference period
nino_off=mean(mean(xa_nino_ens(irm2_a:irm2_b,:),2));

% Annual nino from monthly data (leaves year ends not quite correct)
xa_nino_ann=mon2ann(xa_nino_mn,4,3);
% Annual nino from spatial field
%xa_nino_ann=nino(Xt2m,lat,lon,3,'n');




% Get annual Nino 3.4 value from annual spatial field
efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:33:07.mat';
load([pth,efl],'state_tp','xmeta','Xa_m')
k=find(state_tp=='t');
lat=xmeta{k}.lat;
lon=xmeta{k}.lon;
Xt2m=reshape(Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:),length(lat),length(lon),reconYrs);

% Load mask for Chile/Argentina
%msk=load('ChlArg_mask.mat');
%msk=load('Pampas_mask.mat');
%lndmsk=msk.sftlf;
%lndmsk(lndmsk==0)=NaN;

% Land mask over Chile
msk=load('chile_mask.mat');
lndmsk=msk.chl;


% Latitude weighting
A=cosd(repmat(xmeta{k}.lat,[1 length(xmeta{k}.lon)]));



% Get the ensemble mean; identical to the original reconstruction
Xv_msk=bsxfun(@times,Xt2m,lndmsk);
t2m_casa=wmean_a(Xv_msk,A); % RECONSTRUCTION DROUGHT INDEX




% Drought atlases use JJA instead of annual
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:49:45.mat';
%pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';load([pth,efl])
%pth='/d2/nsteiger/output-da/hydroclimate/'; load([pth,efl])

% PDSI ONLY, ENSEMBLE RECONSTRUCTION
%efl='cesm_lme010_r12000_p998_state_q_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_06-Aug-2018_17:32:52.mat';
%pth='./output-recon/';load([pth,efl],'xmeta')


%=============================================
% PLOT MEAN RECONSTRUCTION WITH UNCERTAINTIES
%=============================================

%load('argentina_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
load('chile_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('chile_megaD_ens_s800.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('chlarg_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
xa_di_sa=Xa_di;
xa_di_ens_sa=Xa_di_ens;
ds_sa=svmd;
mdagr_sa=ens_drght;
sa_yrs=analyrs;

% Reference period
%[~,~,is2] = intersect(analyrs2,sa_yrs); % indices of recon segment

p05 = prctile(xa_di_ens_sa,5,2);
p95 = prctile(xa_di_ens_sa,95,2);


load distinct_colors_20.mat


% PHYDA AND HIGHLIGHT MEGADROUGHTS
figure
hold on
%h1=plot(analyrs3,xa_di_na(ir3),'linewidth',1,'color',[clrs.blue 0.5]); % extra number on the end adds transparency
%h2=plot(analyrs3,xa_di_sa(ir3),'linewidth',1,'color',[clrs.purple 0.5]);
h4=plot(analyrs3,p05(ir3)-mean(xa_di_sa(ir2)),'linewidth',1,'color',clrs.grey);
h5=plot(analyrs3,p95(ir3)-mean(xa_di_sa(ir2)),'linewidth',1,'color',clrs.grey);
h3=plot(analyrs3,xa_di_sa(ir3)-mean(xa_di_sa(ir2)),'linewidth',2,'color',clrs.black);
plot(analyrs3,zeros(length(analyrs3),1),'--','linewidth',2,'color',clrs.yellow)
set(gca,'fontsize',16)
%ylim([-0.8 1])
xlim([analyrs3(1) analyrs3(end)])
xlabel('Year','fontsize',18)
ylabel(['PDSI'],'fontsize',18)
box on
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]

%y=Xa_di(ir3); % using 'l' from above
%yMin=min(y)-0.5;yMax=max(y)+0.5;
ax=gca;yMin=ax.YLim(1);yMax=ax.YLim(2);
ylim([yMin yMax])

% d=ds_na(5:end,:);% only plot droughts that begin after 1000
% % Add a patches for top five times from truth
% for i=1:length(d)
%     p=patch(800+[d(i,1) d(i,2) d(i,2) d(i,1)],[yMin yMin yMax yMax],clrs.blue);
%     p.FaceAlpha=0.35;p.EdgeColor='none';
% end
% 
% d=ds_sa;
% % Add a patches for top five times from truth
% for i=1:length(d)
%     p=patch(analyrs3(1)+[d(i,1) d(i,2) d(i,2) d(i,1)],[yMin yMin yMax yMax],clrs.purple);
%     p.FaceAlpha=0.35;p.EdgeColor='none';
% end

legend([h3,h4],{'PHYDA','5^{th}-95^{th} pctl'},'location','northwest')
hold off

%print(['/home/nsteiger/Dropbox/sa-phyda/drght_phyda_chile_punc.png'],'-dpng','-r300');






%=====================================
% LOAD ENSEMBLE MEGADROUGHT ESTIMATES

% Load both south american and north america drought indices
%load('svmd_ens_90pagree_sa.mat')
%load('svmd_ens_90pagree_pampas.mat')
%load('pampas_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('pampas_megaD_ens_noNAproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
%load('pampas_megaD_ens_noNASWproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')


load('chile_megaD_ens_s800.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
xa_di_sa1=Xa_di;
xa_di_ens_sa1=Xa_di_ens;
ds_sa1=svmd;
mdagr_sa1=ens_drght;

load('chile_megaD_ens_s800_noNAproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
xa_di_sa2=Xa_di;
xa_di_ens_sa2=Xa_di_ens;
ds_sa2=svmd;
mdagr_sa2=ens_drght;

load('chile_megaD_ens_s800_noNASWproxies.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
xa_di_sa3=Xa_di;
xa_di_ens_sa3=Xa_di_ens;
ds_sa3=svmd;
mdagr_sa3=ens_drght;

load('chile_megaD_ens_s800_chileonly.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
xa_di_sa4=Xa_di;
xa_di_ens_sa4=Xa_di_ens;
ds_sa4=svmd;
mdagr_sa4=ens_drght;

%====
% NH
%====

%load('svmd_ens_50pagree.mat')
%load('svmd_ens_90pagree.mat')
load('nasw_megaD_ens.mat','svmd','ens_drght','Xa_di','Xa_di_ens','analyrs','dsens')
xa_di_na=Xa_di;
xa_di_ens_na=Xa_di_ens;
ds_na=svmd;
mdagr_na=ens_drght;

%ds=svmd;




load distinct_colors_20.mat

saoff=mean(xa_di_sa1(ir3));


% PHYDA AND HIGHLIGHT MEGADROUGHTS
figure
hold on
%h1=plot(analyrs3,xa_di_na(ir3),'linewidth',1,'color',[clrs.blue 0.5]); % extra number on the end adds transparency
%h2=plot(analyrs3,xa_di_sa(ir3),'linewidth',1,'color',[clrs.purple 0.5]);
%plot(analyrs3,smooth(DA_di(id3),11),'linewidth',2,'color',clrs.red)
%h3=plot(analyrs3,smooth(xa_di_na(ir3),11),'linewidth',2,'color',clrs.black);
h4=plot(analyrs3,smooth(xa_di_sa1(ir3),11)-saoff,'linewidth',2,'color',clrs.maroon);
h5=plot(analyrs3,smooth(xa_di_sa2(ir3),11)-saoff,'linewidth',2,'color',clrs.navy);
h6=plot(analyrs3,smooth(xa_di_sa3(ir3),11)-saoff,'linewidth',2,'color',clrs.olive);
%h7=plot(analyrs3,smooth(xa_di_sa4(ir3),11),'linewidth',2,'color',clrs.maroon);
plot(analyrs3,zeros(length(analyrs3),1),'--','linewidth',2,'color',clrs.yellow)
set(gca,'fontsize',16)
%ylim([-0.8 1])
xlim([analyrs3(1) analyrs3(end)])
xlabel('Year','fontsize',18)
ylabel(['PDSI'],'fontsize',18)
box on
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 7.5 4.5]; % [left bottom width height]

%y=Xa_di(ir3); % using 'l' from above
%yMin=min(y)-0.5;yMax=max(y)+0.5;
ax=gca;yMin=ax.YLim(1);yMax=ax.YLim(2);
ylim([yMin yMax])

% %----------------------
% % DROUGHT HIGHLIGHTING
% 
% % Droughts in the real data mean
% ds0_n=mdrought(xa_di_na(analyrs),1,'y');
% ds0_s1=mdrought(xa_di_sa1(analyrs),1,'y');
% ds0_s2=mdrought(xa_di_sa2(analyrs),1,'y');
% ds0_s3=mdrought(xa_di_sa3(analyrs),1,'y');
% 
% % Droughts too frequent, throw out based on relative threshold
% %dcut=67;
% dcut=50;
% np67=prctile(ds0_n(:,3),dcut);
% rmv=find(ds0_n(:,3)>np67);
% ds0_n(rmv,:)=[];
% 
% sp67=prctile(ds0_s1(:,3),dcut);
% rmv=find(ds0_s1(:,3)>sp67);
% ds0_s1(rmv,:)=[];
% 
% sp67=prctile(ds0_s2(:,3),dcut);
% rmv=find(ds0_s2(:,3)>sp67);
% ds0_s2(rmv,:)=[];
% 
% sp67=prctile(ds0_s3(:,3),dcut);
% rmv=find(ds0_s3(:,3)>sp67);
% ds0_s3(rmv,:)=[];
% 
% 
% 
% %d=ds_na(5:end,:);% only plot droughts that begin after 1000
% %d=ds_na;
% d=ds0_n;
% for i=1:length(d)
%     p=patch(analyrs(1)+[d(i,1) d(i,2) d(i,2) d(i,1)],[yMin yMin yMax yMax],clrs.blue);
%     p.FaceAlpha=0.35;p.EdgeColor='none';
% end
% 
% %d=ds_sa;
% d=ds0_s1;
% for i=1:length(d)
%     p=patch(analyrs(1)+[d(i,1) d(i,2) d(i,2) d(i,1)],[yMin yMin yMax yMax],clrs.purple);
%     p.FaceAlpha=0.35;p.EdgeColor='none';
% end
% 
% d=ds0_s2;
% for i=1:length(d)
%     p=patch(analyrs(1)+[d(i,1) d(i,2) d(i,2) d(i,1)],[yMin yMin yMax yMax],clrs.orange);
%     p.FaceAlpha=0.35;p.EdgeColor='none';
% end
% 
% d=ds0_s3;
% for i=1:length(d)
%     p=patch(analyrs(1)+[d(i,1) d(i,2) d(i,2) d(i,1)],[yMin yMin yMax yMax],clrs.green);
%     p.FaceAlpha=0.35;p.EdgeColor='none';
% end

%legend([h4,h5,h6,h7],{'Chile','Chile No-NA','Chile No-NASW','Chile Only'},'location','north')
legend([h4,h5,h6],{'SASW','SASW No-NA','SASW No-NASW'},'location','north')
hold off

print(['/home/nsteiger/Dropbox/drght_phyda_chile_11mean_comparisons.png'],'-dpng','-r300');
%print(['/home/nsteiger/Dropbox/sa-phyda/drght_chile_nasw_phyda_11mean_comparisons_wChileOnly.png'],'-dpng','-r300');
%print(['/home/nsteiger/Dropbox/sa-phyda/drght_phyda_chile_11mean_comparisons_highlights.png'],'-dpng','-r300');









