%% BHIz Analyses and Figures

% Plot of time series and scatterplot
year_ins = (1948:1:2010)';
recon_ins_start = find(YEAR == 1948,1,'first');
recon_ins_end = find(YEAR == 2009,1,'last');
panel = {'a','b','c','d'};%,'D'};
figure
subplot(3,2,2)%2,2,[1,2])
ax = gca;
waterblue = rgb('water blue');
plot(year_ins(1:62),BHIz_ins(1:62),'k-','LineWidth',3)
hold on
plot(year_ins(1:62),RECON(recon_ins_start:recon_ins_end),'-','Color',waterblue,'LineWidth',3)
hold off
hold on
patch([year_ins(1:62); flipud(year_ins(1:62))],[PREUP_95(recon_ins_start:recon_ins_end); flipud(RECON(recon_ins_start:recon_ins_end))],waterblue,'FaceAlpha',0.1,'LineStyle','none')
patch([year_ins(1:62); flipud(year_ins(1:62))],[PRELO_05(recon_ins_start:recon_ins_end); flipud(RECON(recon_ins_start:recon_ins_end))],waterblue,'FaceAlpha',0.1,'LineStyle','none')
hold off
axis tight
ax.YColor = 'k';
ax.XColor = 'k';
ax.FontSize = 20;
ax.FontWeight = 'bold';
ax.XMinorTick = 'on'; ax.YMinorTick = 'on';
ax.XGrid = 'on'; ax.YGrid = 'on';
ylabel('BHIz','FontSize',20,'FontWeight','bold')
xlim([1948 2010])
xlabel('Year (CE)','FontSize',20,'FontWeight','bold')
legend('Instrumental BHI','Reconstructed BHI')
ntitle(panel{2},'location','nw','fontsize',20,'fontweight','bold')
subplot(3,2,1)%2,2,[3,4])
ax = gca;
greyblue = rgb('grey/blue');
scatter(BHIz_ins(1:62),RECON(recon_ins_start:recon_ins_end),56,greyblue,'filled')
hold on
ylim([-2 2])
plot([ylim ylim],[ylim ylim],'k--','LineWidth',2)
hold off
% hline = refline;
% hline.LineStyle = '--';
% hline.Color = 'k';
% hline.LineWidth = 2;
ax.Box = 'on';
ax.XMinorTick = 'on'; ax.YMinorTick = 'on';
ax.XGrid = 'on'; ax.YGrid = 'on';
ax.FontSize = 20;
ax.FontWeight = 'bold';
ylabel('Reconstructed BHIz','FontSize',20,'FontWeight','bold')
xlabel('Instrumental BHIz','FontSize',20,'FontWeight','bold')
ntitle(panel{1},'location','nw','fontsize',20,'fontweight','bold')
subplot(3,2,[3,4])
% ax = gca;
% waterblue = rgb('water blue');
% plot(YEAR,RECON,'-','Color',waterblue,'LineWidth',1.5)
% % hold on
% % plot(YEAR,BHIz_recon_10movmean,'-','Color',midnight,'LineWidth',3)
% % hold off
% hold on
% patch([YEAR; flipud(YEAR)],[PREUP_95; flipud(PRELO_05)],waterblue,'FaceAlpha',0.3,'LineStyle','none')
% hold off
% axis tight
% ax.YColor = 'k';
% ax.XColor = 'k';
% ylabel('Reconstructed BHI')
% xlim([1430 2011])
% ntitle(panel{3},'location','nw')
%subplot(3,2,[5,6])
ax = gca;
forest = rgb('forest');
%burntsienna = rgb('burnt sienna');
yyaxis left
plot(YEAR,CORES,'-','Color',waterblue,'LineWidth',2)
ax.YColor = 'k';%burntsienna;
ax.YMinorTick = 'on';
ax.FontSize = 20;
ax.FontWeight = 'bold';
ylabel({'Number of ';'Cores'},'FontSize',20,'FontWeight','bold')
axis tight
yyaxis right
plot(YEAR,CRNS,'-','Color',forest,'LineWidth',2)
ax.YColor = 'k';%forest;
ax.XMinorTick = 'on'; ax.YMinorTick = 'on';
ylabel({'Number of ';'Chronologies'},'FontSize',20,'FontWeight','bold','Rotation',270)
ax.FontSize = 20;
ax.FontWeight = 'bold';
axis tight
%xlabel('Year (CE)')
xlim([1140 2010])
legend('Number of Cores','Number of Chronologies','Location','north')
ntitle(panel{3},'location','nw','fontsize',20,'fontweight','bold')
subplot(3,2,[5,6])
ax = gca;
yyaxis left
plot(YEAR,CVRE,'k-','LineWidth',2)
hold on
plot(YEAR,CRSQ,'r-','LineWidth',2)
hold off
ax.YColor = 'k';
ax.YMinorTick = 'on';
ax.FontSize = 20;
ax.FontWeight = 'bold';
ylim([0 0.8]);
ax.YTick = 0.1:0.1:0.9;
ylabel('CVRE/CRSQ','FontSize',20,'FontWeight','bold')
axis tight
yyaxis right
%plot(YEAR,CRSQ,'r-','LineWidth',2)
plot(YEAR,PCs,'-','Color',waterblue,'LineWidth',2)
ax.YColor = 'k';
ax.YMinorTick = 'on';
ax.XMinorTick = 'on';
ylim([0 12]); %ylim([0 0.8]);
ax.YTick = 0:2:12;%ax.YTick = 0.2:0.1:0.7;
ax.FontSize = 20;
ax.FontWeight = 'bold';
ylabel({'Number'; 'of PCs'},'FontSize',20,'FontWeight','bold','Rotation',270);%ylabel('CRSQ','FontSize',20,'FontWeight','bold','Rotation',270)
xlabel('Year (CE)','FontSize',20,'FontWeight','bold')
legend('CVRE','CRSQ','Number of PCs','Location','nw')
ntitle(panel{4},'location','nw','fontsize',20,'fontweight','bold')
% subplot(4,2,[7,8])
% ax = gca;
% yyaxis left
% plot(YEAR,PCs,'k-','LineWidth',2)
% ax.YColor = 'k';
% ax.YMinorTick = 'on';
% ax.FontSize = 20;
% ax.FontWeight = 'bold';
% ylim([0 12]);
% ax.YTick = 0:2:12;
% ylabel({'Number';'of PCs'},'FontSize',20,'FontWeight','bold')
% axis tight
% yyaxis right
% plot(YEAR,Nests,'r-','LineWidth',2)
% ax.YColor = 'k';
% ax.YMinorTick = 'on';
% ax.XMinorTick = 'on';
% ylim([0 50]);
% ax.YTick = 0:10:50;
% ax.FontSize = 20;
% ax.FontWeight = 'bold';
% ylabel('Nests','FontSize',20,'FontWeight','bold','Rotation',270)
% xlabel('Year (CE)','FontSize',20,'FontWeight','bold')
% legend('Number of PCs','Nests','Location','nw')
% ntitle(panel{5},'location','nw','fontsize',20,'fontweight','bold')

%
%% Validation statistics

% Mean absolute error
bhi_mae = sum((RECON(find(YEAR == 1948,1,'first'):find(YEAR == 2009,1,'first'))-BHIz_ins(find(year_ins == 1948,1,'first'):find(year_ins == 2009,1,'first'))))/length(1948:2009)

%
%% F-tests

[mwp20_h, mwp20_p, mwp20_ci, mwp20_stats] = vartest2(RECON(find(YEAR == 1900,1,'first'):...
    find(YEAR == 2009,1,'first')),RECON(1:find(YEAR == 1250,1,'first')),'Tail','right');

[mwpind_h, mwpind_p, mwpind_ci, mwpind_stats] = vartest2(RECON(find(YEAR == 1851,1,'first'):...
    find(YEAR == 2009,1,'first')),RECON(1:find(YEAR == 1250,1,'first')),'Tail','right');

[lia20_h, lia20_p, lia20_ci, lia20_stats] = vartest2(RECON(find(YEAR == 1900,1,'first'):...
    find(YEAR == 2009,1,'first')),RECON(find(YEAR == 1300,1,'first'):find(YEAR == 1850,1,'first')),...
    'Tail','right');

[pre20_h, pre20_p, pre20_ci, pre20_stats] = vartest2(RECON(find(YEAR == 1900,1,'first'):...
    find(YEAR == 2009,1,'first')),RECON(1:find(YEAR == 1850,1,'first')),'Tail','right')

%

BHIz_recon_10movmean = movmean(RECON,10,'Endpoints','fill');
BHIz_recon_pup95_10movmean = movmean(PREUP_95,10,'Endpoints','fill');
BHIz_recon_plo05_10movmean = movmean(PRELO_05,10,'Endpoints','fill');
BHIz_recon_5movmean = movmean(RECON,5,'Endpoints','fill');
BHIz_recon_pup95_5movmean = movmean(PREUP_95,5,'Endpoints','fill');
BHIz_recon_plo05_5movmean = movmean(PRELO_05,5,'Endpoints','fill');
BHIz_recon_30movmean = movmean(RECON,30,'Endpoints','fill');
BHIz_recon_pup95_30movmean = movmean(PREUP_95,30,'Endpoints','fill');
BHIz_recon_plo05_30movmean = movmean(PRELO_05,30,'Endpoints','fill');

%

panel = {'a','b','c','d','e','f'};%,'g','h'};
%midnight = rgb('midnight');
%ultramar = rgb('ultramarine');
BHIz_ins_10mm = movmean(BHIz_ins,10,'EndPoints','fill');
%BHIz_ins_5mm = movmean(BHIz_ins,5,'EndPoints','fill');
BHIz_ins_30mm = movmean(BHIz_ins,30,'Endpoints','fill');
figure
%subplot(6,2,[1,2])
ax(1) = subsubplot(6,1,1,'vpad',0.01);
%ax = gca;
waterblue = rgb('water blue');
plot(YEAR,RECON,'-','Color',waterblue,'LineWidth',2)
% hold on
% plot(YEAR,BHIz_recon_10movmean,'-','Color',midnight,'LineWidth',3)
% hold off
hold on
patch([YEAR; flipud(YEAR)],[PREUP_95; flipud(PRELO_05)],waterblue,'FaceAlpha',0.3,'LineStyle','none')
hold off
axis tight
ylabel('BHI','FontSize',20,'FontWeight','bold')
ax(1).XAxisLocation = 'top';
ax(1).XColor = 'k';
ax(1).FontWeight = 'bold';
ax(1).FontSize = 15;
ax(1).YMinorTick = 'on';
ax(1).XMinorTick = 'on';
ax(1).XTick = 1150:50:2000;
ax(1).XTickLabel = {'','1200','','1300','','1400','','1500','','1600','','1700','','1800','','1900','','2000'};
ntitle(panel{1},'location','nw','fontsize',13,'fontweight','bold')
%subplot(6,2,[9,10])
% Originally axis 5
ax(3) = subsubplot(6,1,3,'vpad',0.01);
%ax = gca;
plot(YEAR,BHIz_recon_10movmean,'-','Color',waterblue,'LineWidth',3)
hold on
patch([YEAR(find(BHIz_recon_10movmean > -100,1,'first'):find(BHIz_recon_10movmean > -100,1,'last'));...
    flipud(YEAR(find(BHIz_recon_10movmean > -100,1,'first'):find(BHIz_recon_10movmean > -100,1,'last')))],...
    [BHIz_recon_pup95_10movmean(find(BHIz_recon_10movmean > -100,1,'first'):find(BHIz_recon_10movmean > -100,1,'last'));...
    flipud(BHIz_recon_plo05_10movmean(find(BHIz_recon_10movmean > -100,1,'first'):find(BHIz_recon_10movmean > -100,1,'last')))],...
    waterblue,'FaceAlpha',0.3,'LineStyle','none')
hold off
axis tight
ylabel('BHI (10-yr mean)','FontSize',20,'FontWeight','bold')
ax(3).YMinorTick = 'on';
ax(3).FontSize = 15;
ax(3).FontWeight = 'bold';
xlim([1140 2009])
ntitle(panel{3},'location','nw','fontsize',13,'fontweight','bold')
%subplot(6,2,[3,4])
ax(2) = subsubplot(6,1,2,'vpad',0.01);
%anomaly(YEAR,RECON,'thresh',mean(BHIz_ins,'omitnan'))%BHIz_recon_10movmean)
anomaly(YEAR,RECON,'base',mean(BHIz_ins,'omitnan'))
axis tight
%subplot(6,2,[11,12])
ylabel('BHI Anomaly','FontSize',20,'FontWeight','bold','Rotation',270)
ax(2).FontWeight = 'bold';
ax(2).FontSize = 15;
ax(2).YMinorTick = 'on';
ntitle(panel{2},'location','nw','fontsize',13,'fontweight','bold')
% Originally axis 6
ax(4) = subsubplot(6,1,4,'vpad',0.01)
%anomaly(YEAR,BHIz_recon_10movmean,'thresh',mean(BHIz_ins_10mm,'omitnan'))
anomaly(YEAR,BHIz_recon_10movmean,'base',mean(BHIz_ins,'omitnan'))%mean(BHIz_ins_10mm,'omitnan'))
axis tight
ylabel({'BHI Anomaly'; '(10-yr mean)'},'FontSize',20,'FontWeight','bold','Rotation',270)
ax(4).FontWeight = 'bold';
ax(4).FontSize = 15;
ax(4).YMinorTick = 'on';
xlim([1140 2009])
% ax(4).XMinorTick = 'on';
% ax(4).XTick = 1150:50:2000;
% xlabel('Year (CE)','FontSize',15,'FontWeight','bold')
ntitle(panel{4},'location','nw','fontsize',13,'fontweight','bold')
%subplot(6,2,[5,6])
% Originally axis 3 and 5-yr mm
ax(5) = subsubplot(6,1,5,'vpad',0.01);
%ax = gca;
plot(YEAR,BHIz_recon_30movmean,'-','Color',waterblue,'LineWidth',3)
hold on
patch([YEAR(find(BHIz_recon_30movmean > -100,1,'first'):find(BHIz_recon_30movmean > -100,1,'last'));...
    flipud(YEAR(find(BHIz_recon_30movmean > -100,1,'first'):find(BHIz_recon_30movmean > -100,1,'last')))],...
    [BHIz_recon_pup95_30movmean(find(BHIz_recon_30movmean > -100,1,'first'):find(BHIz_recon_30movmean > -100,1,'last'));...
    flipud(BHIz_recon_plo05_30movmean(find(BHIz_recon_30movmean > -100,1,'first'):find(BHIz_recon_30movmean > -100,1,'last')))],...
    waterblue,'FaceAlpha',0.3,'LineStyle','none')
hold off
axis tight
ylabel('BHI (30-yr mean)','FontSize',20,'FontWeight','bold')
ax(5).YMinorTick = 'on';
ax(5).FontWeight = 'bold';
ax(5).FontSize = 15;
xlim([1140 2009])
ntitle(panel{5},'location','nw','fontsize',13,'fontweight','bold')
%subplot(6,2,[7,8])
% Originally axis 4 and 5-yr mm
ax(6) = subsubplot(6,1,6,'vpad',0.01);
%ax = gca;
%anomaly(YEAR,BHIz_recon_30movmean,'thresh',mean(BHIz_ins_30mm,'omitnan'))
anomaly(YEAR,BHIz_recon_30movmean,'base',mean(BHIz_ins,'omitnan'))%mean(BHIz_ins_30mm,'omitnan'))
axis tight
ylabel({'BHI Anomaly'; '(30-yr mean)'},'FontSize',20,'FontWeight','bold','Rotation',270)
%xlabel('Year (CE)','FontSize',15,'FontWeight','bold')
ax(6).YMinorTick = 'on';
ax(6).FontSize = 15;
ax(6).FontWeight = 'bold';
xlim([1140 2009])
ax(6).XMinorTick = 'on';
ax(6).XTick = 1150:50:2000;
ax(6).XTickLabel = {'','1200','','1300','','1400','','1500','','1600','','1700','','1800','','1900','','2000'};
xlabel('Year (CE)','FontSize',20,'FontWeight','bold')
ntitle(panel{6},'location','nw','fontsize',15,'fontweight','bold')
% ax(5) = subsubplot(8,1,5,'vpad',0.01);
% plot(YEAR,recon_5mv,'k-','LineWidth',3)
% axis tight
% ylabel({'Variance'; '(5-yr)'},'FontSize',15,'FontWeight','bold')
% ax(5).YMinorTick = 'on';
% ax(5).FontWeight = 'bold';
% ax(5).FontSize = 13;
% ntitle(panel{5},'location','nw','fontsize',13,'fontweight','bold')
% ax(8) = subsubplot(8,1,8,'vpad',0.01);
% plot(YEAR,recon_10mv,'k-','LineWidth',3)
% axis tight
% ylabel({'Variance'; '(10-yr)'},'FontSize',15,'FontWeight','bold','Rotation',270)
% ax(8).FontWeight = 'bold';
% ax(8).FontSize = 13;
% ax(8).YMinorTick = 'on';
% ax(8).XMinorTick = 'on';
% ax(8).XTick = 1150:50:2000;
% xlabel('Year (CE)','FontSize',15,'FontWeight','bold')
% ntitle(panel{8},'location','nw','fontsize',13,'fontweight','bold')
set(ax(2:2:6),'YAxisLocation','right')

%
%% Testing out some other moving means

BHIz_recon_15movmean = movmean(RECON,15,'Endpoints','fill');
BHIz_recon_pup95_15movmean = movmean(PREUP_95,15,'Endpoints','fill');
BHIz_recon_plo05_15movmean = movmean(PRELO_05,15,'Endpoints','fill');
BHIz_ins_15mm = movmean(BHIz_ins,15,'EndPoints','fill');

BHIz_recon_40movmean = movmean(RECON,40,'Endpoints','fill');
BHIz_recon_pup95_40movmean = movmean(PREUP_95,40,'Endpoints','fill');
BHIz_recon_plo05_40movmean = movmean(PRELO_05,40,'Endpoints','fill');
BHIz_ins_40mm = movmean(BHIz_ins,40,'Endpoints','fill');

BHIz_recon_25movmean = movmean(RECON,25,'Endpoints','fill');
BHIz_recon_pup95_25movmean = movmean(PREUP_95,25,'Endpoints','fill');
BHIz_recon_plo05_25movmean = movmean(PRELO_05,25,'Endpoints','fill');
BHIz_ins_25mm = movmean(BHIz_ins,25,'Endpoints','fill');

%

% Plotting the moving means from the section above.
panel = {'a','b','c','d','e','f'};
figure
%subplot(6,2,[1,2])
ax(5) = subsubplot(6,1,5,'vpad',0.01);
%ax = gca;
waterblue = rgb('water blue');
plot(YEAR,BHIz_recon_40movmean,'-','Color',waterblue,'LineWidth',2)
% hold on
% plot(YEAR,BHIz_recon_10movmean,'-','Color',midnight,'LineWidth',3)
% hold off
hold on
patch([YEAR(find(BHIz_recon_40movmean > -100, 1, 'first'):find(BHIz_recon_40movmean > -100,1,'last'));...
    flipud(YEAR(find(BHIz_recon_40movmean > -100, 1, 'first'):find(BHIz_recon_40movmean > -100,1,'last')))],...
    [BHIz_recon_pup95_40movmean(find(BHIz_recon_40movmean > -100, 1, 'first'):find(BHIz_recon_40movmean > -100,1,'last'));...
    flipud(BHIz_recon_plo05_40movmean(find(BHIz_recon_40movmean > -100, 1, 'first'):find(BHIz_recon_40movmean > -100,1,'last')))],...
    waterblue,'FaceAlpha',0.3,'LineStyle','none')
hold off
axis tight
xlim([1140 2010]);
ylabel('BHI (40-yr mean)','FontSize',15,'FontWeight','bold')
ax(5).YMinorTick = 'on';
ax(5).FontSize = 13;
ax(5).FontWeight = 'bold';
% ax(1).XAxisLocation = 'top';
% ax(1).XColor = 'k';
ntitle(panel{5},'location','nw','fontsize',10,'fontweight','bold')
%subplot(6,2,[9,10])

ax(1) = subsubplot(6,1,1,'vpad',0.01);
%ax = gca;
plot(YEAR,BHIz_recon_15movmean,'-','Color',waterblue,'LineWidth',3)
hold on
patch([YEAR(find(BHIz_recon_15movmean > -100,1,'first'):find(BHIz_recon_15movmean > -100,1,'last'));...
    flipud(YEAR(find(BHIz_recon_15movmean > -100,1,'first'):find(BHIz_recon_15movmean > -100,1,'last')))],...
    [BHIz_recon_pup95_15movmean(find(BHIz_recon_15movmean > -100,1,'first'):find(BHIz_recon_15movmean > -100,1,'last'));...
    flipud(BHIz_recon_plo05_15movmean(find(BHIz_recon_15movmean > -100,1,'first'):find(BHIz_recon_15movmean > -100,1,'last')))],...
    waterblue,'FaceAlpha',0.3,'LineStyle','none')
hold off
axis tight
ylabel('BHI (15-yr mean)','FontSize',15,'FontWeight','bold')
ax(1).XAxisLocation = 'top';
ax(1).XColor = 'k';
ax(1).YMinorTick = 'on';
ax(1).XMinorTick = 'on';
ax(1).FontSize = 13;
ax(1).FontWeight = 'bold';
ax(1).XTick = 1150:50:2000;
ntitle(panel{1},'location','nw','fontsize',10,'fontweight','bold')
%subplot(6,2,[3,4])

ax(6) = subsubplot(6,1,6,'vpad',0.01);
anomaly(YEAR,BHIz_recon_40movmean,'thresh',mean(BHIz_ins_40mm,'omitnan'))%BHIz_recon_10movmean)
axis tight
%subplot(6,2,[11,12])
ylabel(['BHI Anomaly '; '(40-yr mean)'],'Fontsize',15,'FontWeight','bold','Rotation',270)
xlim([1140 2009])
xlabel('Year (CE)','FontSize',15,'FontWeight','bold')
ax(6).YMinorTick = 'on';
ax(6).XMinorTick = 'on';
ax(6).XTick = 1150:50:2000;
ax(6).FontSize = 13;
ax(6).FontWeight = 'bold';
ntitle(panel{6},'location','nw','fontsize',10,'fontweight','bold')

ax(2) = subsubplot(6,1,2,'vpad',0.01)
anomaly(YEAR,BHIz_recon_15movmean,'thresh',mean(BHIz_ins_15mm,'omitnan'))
axis tight
ylabel(['BHI Anomaly ';'(15-yr mean)'],'FontSize',15,'FontWeight','bold','Rotation',270)
ax(2).YMinorTick = 'on';
ax(2).FontSize = 13;
ax(2).FontWeight = 'bold';
ntitle(panel{2},'location','nw','fontsize',10,'fontweight','bold')
%subplot(6,2,[5,6])

ax(3) = subsubplot(6,1,3,'vpad',0.01);
%ax = gca;
plot(YEAR,BHIz_recon_25movmean,'-','Color',waterblue,'LineWidth',3)
hold on
patch([YEAR(find(BHIz_recon_25movmean > -100,1,'first'):find(BHIz_recon_25movmean > -100,1,'last'));...
    flipud(YEAR(find(BHIz_recon_25movmean > -100,1,'first'):find(BHIz_recon_25movmean > -100,1,'last')))],...
    [BHIz_recon_pup95_25movmean(find(BHIz_recon_25movmean > -100,1,'first'):find(BHIz_recon_25movmean > -100,1,'last'));...
    flipud(BHIz_recon_plo05_25movmean(find(BHIz_recon_25movmean > -100,1,'first'):find(BHIz_recon_25movmean > -100,1,'last')))],...
    waterblue,'FaceAlpha',0.3,'LineStyle','none')
hold off
axis tight
ylabel('BHI (25-yr mean)','FontSize',15,'FontWeight','bold')
ax(3).YMinorTick = 'on';
ax(3).FontSize = 13;
ax(3).FontWeight = 'bold';
ntitle(panel{3},'location','nw','fontsize',10,'fontweight','bold')
%subplot(6,2,[7,8])

ax(4) = subsubplot(6,1,4,'vpad',0.01);
%ax = gca;
anomaly(YEAR,BHIz_recon_25movmean,'thresh',mean(BHIz_ins_25mm,'omitnan'))
axis tight
ylabel(['BHI Anomaly ';'(25-yr mean)'],'FontSize',15,'FontWeight','bold','Rotation',270)
ax(4).YMinorTick = 'on';
ax(4).FontSize = 13;
ax(4).FontWeight = 'bold';
%xlabel('Year (CE)')
ntitle(panel{4},'location','nw','fontsize',10,'fontweight','bold')
set(ax(2:2:6),'YAxisLocation','right')

%
%% Instrumental Data Extraction

% pdsi_cru = ncread('scPDSI.cru_ts4.05early1.1901.2020.cal_1901_20.bams.2021.GLOBAL.1901.2020.nc','scpdsi');
% pdsi_time_cru = ncread('scPDSI.cru_ts4.05early1.1901.2020.cal_1901_20.bams.2021.GLOBAL.1901.2020.nc','time');
% pdsi_lat_cru = ncread('scPDSI.cru_ts4.05early1.1901.2020.cal_1901_20.bams.2021.GLOBAL.1901.2020.nc','latitude');
% pdsi_lon_cru = ncread('scPDSI.cru_ts4.05early1.1901.2020.cal_1901_20.bams.2021.GLOBAL.1901.2020.nc','longitude');
% 
% pdsi_newtime_cru = datetime(datestr(daysadd('01-Jan-1900',pdsi_time_cru)));
% Mayst = find(pdsi_newtime_cru == '01-May-1948',1,'first');
% Junst = find(pdsi_newtime_cru == '01-Jun-1948',1,'first');
% Julst = find(pdsi_newtime_cru == '01-Jul-1948',1,'first');
% lated = find(pdsi_lat_cru >= 26,1,'last');
% latst = find(pdsi_lat_cru <= 38,1,'first');
% lonst = find(pdsi_lon_cru >= -92,1,'first');
% loned = find(pdsi_lon_cru <= -70,1,'last');
% regional_MJJ_pdsi_cru = nan(length(lonst:loned),length(latst:lated),length(1948:2010),3);
% regional_pdsi_cru = pdsi_cru(lonst:loned,latst:lated,:);
% 
% for i = 1:length(1948:2010)
%     regional_MJJ_pdsi_cru(:,:,i,1) = regional_pdsi_cru(:,:,Mayst);
%     regional_MJJ_pdsi_cru(:,:,i,2) = regional_pdsi_cru(:,:,Junst);
%     regional_MJJ_pdsi_cru(:,:,i,3) = regional_pdsi_cru(:,:,Julst);
%     Mayst = Mayst+12;
%     Junst = Junst+12;
%     Julst = Julst+12;
% end
% 
% pdsi_monmean = ncread('pdsi.mon.mean.selfcalibrated.nc','pdsi');
% pdsi_lat_monmean = ncread('pdsi.mon.mean.selfcalibrated.nc','lat');
% pdsi_lon_monmean = ncread('pdsi.mon.mean.selfcalibrated.nc','lon');
% pdsi_time_monmean = ncread('pdsi.mon.mean.selfcalibrated.nc','time');
% 
% pdsi_monmean_newtime = datetime(datestr(datenum(1800,1,1,double(pdsi_time_monmean),0,0)));
% Mayst = find(pdsi_monmean_newtime == '01-May-1948',1,'first');
% Junst = find(pdsi_monmean_newtime == '01-Jun-1948',1,'first');
% Julst = find(pdsi_monmean_newtime == '01-Jul-1948',1,'first');
% latst = find(pdsi_lat_monmean >= 26,1,'first');
% lated = find(pdsi_lat_monmean <= 38,1,'last');
% lonst = find(pdsi_lon_monmean >= -92,1,'first');
% loned = find(pdsi_lon_monmean <= -70,1,'last');
% regional_pdsi_monmean = pdsi_monmean(lonst:loned,latst:lated,:);
% regional_MJJ_pdsi_monmean = nan(length(lonst:loned),length(latst:lated),length(1948:2010),3);
% 
% for i = 1:length(1948:2010)
%     regional_MJJ_pdsi_monmean(:,:,i,1) = regional_pdsi_monmean(:,:,Mayst);
%     regional_MJJ_pdsi_monmean(:,:,i,2) = regional_pdsi_monmean(:,:,Junst);
%     regional_MJJ_pdsi_monmean(:,:,i,3) = regional_pdsi_monmean(:,:,Julst);
%     Mayst = Mayst+12;
%     Junst = Junst+12;
%     Julst = Julst+12;
% end
% 
% precip_monmean_1x1 = ncread('precip.mon.mean.1x1.nc','precip');
% precip_monmean_1x1_time = ncread('precip.mon.mean.1x1.nc','time');
% precip_monmean_1x1_lon = ncread('precip.mon.mean.1x1.nc','lon');
% precip_monmean_1x1_lat = ncread('precip.mon.mean.1x1.nc','lat');
% 
% precip_monmean_1x1_newtime = datetime(datestr(datenum(1800,1,1,double(precip_monmean_1x1_time),0,0)));
% Mayst = find(precip_monmean_1x1_newtime == '01-May-1948',1,'first');
% Junst = find(precip_monmean_1x1_newtime == '01-Jun-1948',1,'first');
% Julst = find(precip_monmean_1x1_newtime == '01-Jul-1948',1,'first');
% latst = find(precip_monmean_1x1_lat <= 38,1,'first');
% lated = find(precip_monmean_1x1_lat >= 26,1,'last');
% lonst = find((precip_monmean_1x1_lon-180) >= 70,1,'first');
% loned = find((precip_monmean_1x1_lon-180) <= 92,1,'last');
% regional_precip_monmean_1x1 = precip_monmean_1x1(lonst:loned,latst:lated,:);
% regional_MJJ_precip_monmean_1x1 = nan(length(lonst:loned),length(latst:lated),length(1948:2010),3);
% 
% for i = 1:length(1948:2010)
%     regional_MJJ_precip_monmean_1x1(:,:,i,1) = regional_precip_monmean_1x1(:,:,Mayst);
%     regional_MJJ_precip_monmean_1x1(:,:,i,2) = regional_precip_monmean_1x1(:,:,Junst);
%     regional_MJJ_precip_monmean_1x1(:,:,i,3) = regional_precip_monmean_1x1(:,:,Julst);
%     Mayst = Mayst+12;
%     Junst = Junst+12;
%     Julst = Julst+12;
% end
% 
% soilw_monmean = ncread('soilw.mon.mean.v2.nc','soilw');
% soilw_monmean_time = ncread('soilw.mon.mean.v2.nc','time');
% soilw_monmean_lat = ncread('soilw.mon.mean.v2.nc','lat');
% soilw_monmean_lon = ncread('soilw.mon.mean.v2.nc','lon');
% 
% soilw_monmean_newtime = datetime(datestr(datenum(1800,1,double(soilw_monmean_time),0,0,0)));
% Mayst = find(soilw_monmean_newtime == '31-May-1948',1,'first');
% Junst = find(soilw_monmean_newtime == '30-Jun-1948',1,'first');
% Julst = find(soilw_monmean_newtime == '31-Jul-1948',1,'first');
% latst = find(soilw_monmean_lat <= 38,1,'first');
% lated = find(soilw_monmean_lat >= 26,1,'last');
% lonst = find((soilw_monmean_lon-180) >= 70,1,'first');
% loned = find((soilw_monmean_lon-180) <= 92,1,'last');
% regional_soilw_monmean = soilw_monmean(lonst:loned,latst:lated,:);
% regional_MJJ_soilw_monmean = nan(length(lonst:loned),length(latst:lated),length(1948:2010),3);
% 
% for i = 1:length(1948:2010)
%     regional_MJJ_soilw_monmean(:,:,i,1) = regional_soilw_monmean(:,:,Mayst);
%     regional_MJJ_soilw_monmean(:,:,i,2) = regional_soilw_monmean(:,:,Junst);
%     regional_MJJ_soilw_monmean(:,:,i,3) = regional_soilw_monmean(:,:,Julst);
%     Mayst = Mayst+12;
%     Junst = Junst+12;
%     Julst = Julst+12;
% end
% 
% precip_montot_1x1 = ncread('precip.mon.total.1x1.v2018.nc','precip');
% precip_montot_1x1_time = ncread('precip.mon.total.1x1.v2018.nc','time');
% precip_montot_1x1_lon = ncread('precip.mon.total.1x1.v2018.nc','lon');
% precip_montot_1x1_lat = ncread('precip.mon.total.1x1.v2018.nc','lat');
% 
% precip_montot_1x1_newtime = datetime(datestr(datenum(1800,1,double(precip_montot_1x1_time),0,0,0)));
% Mayst = find(precip_montot_1x1_newtime == '31-May-1948',1,'first');
% Junst = find(precip_montot_1x1_newtime == '30-Jun-1948',1,'first');
% Julst = find(precip_montot_1x1_newtime == '31-Jul-1948',1,'first');
% latst = find(precip_montot_1x1_lat <= 38,1,'first');
% lated = find(precip_montot_1x1_lat >= 26,1,'last');
% lonst = find((precip_montot_1x1_lon-180) >= 70,1,'first');
% loned = find((precip_montot_1x1_lon-180) <= 92,1,'last');
% regional_precip_montot_1x1 = precip_montot_1x1(lonst:loned,latst:lated,:);
% regional_MJJ_precip_montot_1x1 = nan(length(lonst:loned),length(latst:lated),length(1948:2010),3);
% 
% for i = 1:length(1948:2010)
%     regional_MJJ_precip_montot_1x1(:,:,i,1) = regional_precip_montot_1x1(:,:,Mayst);
%     regional_MJJ_precip_montot_1x1(:,:,i,2) = regional_precip_montot_1x1(:,:,Junst);
%     regional_MJJ_precip_montot_1x1(:,:,i,3) = regional_precip_montot_1x1(:,:,Julst);
%     Mayst = Mayst+12;
%     Junst = Junst+12;
%     Julst = Julst+12;
% end

%
%% Instrumental Data Comparisons

% This was first edited on July 14, 2022

% regional_MJJ_pdsi_cru(regional_MJJ_pdsi_cru < -100) = NaN;
% regional_MJJ_pdsi_cru_gridded_mean = squeeze(mean(regional_MJJ_pdsi_cru,4,'omitnan'));
% lated = find(pdsi_lat_cru >= 26,1,'last');
% latst = find(pdsi_lat_cru <= 38,1,'first');
% lonst = find(pdsi_lon_cru >= -92,1,'first');
% loned = find(pdsi_lon_cru <= -70,1,'last');
% regional_MJJ_pdsi_cru_lat = pdsi_lat_cru(latst:lated);
% regional_MJJ_pdsi_cru_lon = pdsi_lon_cru(lonst:loned);
% regional_MJJ_pdsi_cru_latmat = repmat(regional_MJJ_pdsi_cru_lat,[1 length(regional_MJJ_pdsi_cru_lon)])';
% regional_MJJ_pdsi_cru_weighted_mean = nan(1,length(1948:2010))';%wmean(regional_MJJ_cru_grided_mean,cosd(regional_MJJ_pdsi_cru_latmat),'omitnan'));
% for i = 1:length(1948:2010)
%     regional_MJJ_pdsi_cru_weighted_mean(i,1) = wmean(...
%         squeeze(regional_MJJ_pdsi_cru_gridded_mean(:,:,i)),...
%         cosd(regional_MJJ_pdsi_cru_latmat),'all','omitnan');
% end
% 
% regional_MJJ_pdsi_cru_weighted_mean_detrended = detrend(regional_MJJ_pdsi_cru_weighted_mean);
% % Spatially-weighted correlations
% [rins_pdsi_cru, pins_pdsi_cru] = corr(detrend(BHIz_ins),regional_MJJ_pdsi_cru_weighted_mean_detrended,...
%     'Type','Pearson','Rows','pairwise','Tail','right');
% [rrecon_pdsi_cru, precon_pdsi_cru] = corr(detrend(RECON(find(YEAR == 1948,1,'first'):find(YEAR == 2009,1,'first'))),...
%     regional_MJJ_pdsi_cru_weighted_mean_detrended(1:62),'Type','Pearson','Rows','pairwise','Tail','right');
% % Gridded correlations
% [gx,gy,~] = size(regional_MJJ_pdsi_cru_gridded_mean);
% rins_pdsi_cru_gridded = nan([gx gy]);
% pins_pdsi_cru_gridded = nan([gx gy]);
% rrecon_pdsi_cru_gridded = nan([gx gy]);
% precon_pdsi_cru_gridded = nan([gx gy]);
% for i = 1:gx
%     for j = 1:gy
%         pdsi_cru_ts = squeeze(regional_MJJ_pdsi_cru_gridded_mean(i,j,:));
%         [r,p] = corr(detrend(BHIz_ins),detrend(pdsi_cru_ts),'Type','Pearson','Rows','pairwise','Tail','right');
%         rins_pdsi_cru_gridded(i,j) = r;
%         pins_pdsi_cru_gridded(i,j) = p;
%         [r_r,p_r] = corr(detrend(RECON(find(YEAR == 1948,1,'first'):find(YEAR == 2009,1,'first'))),...
%             detrend(pdsi_cru_ts(1:62)),'Type','Pearson','Rows','pairwise','Tail','right');
%         rrecon_pdsi_cru_gridded(i,j) = r_r;
%         precon_pdsi_cru_gridded(i,j) = p_r;
%     end
% end
% % Start making correlation maps here
% 
% regional_MJJ_pdsi_monmean(regional_MJJ_pdsi_monmean < -100) = NaN;
% regional_MJJ_pdsi_monmean_gridded_mean = squeeze(mean(regional_MJJ_pdsi_monmean,4,'omitnan'));
% latst = find(pdsi_lat_monmean >= 26,1,'first');
% lated = find(pdsi_lat_monmean <= 38,1,'last');
% lonst = find(pdsi_lon_monmean >= -92,1,'first');
% loned = find(pdsi_lon_monmean <= -70,1,'last');
% regional_MJJ_pdsi_monmean_lat = pdsi_lat_monmean(latst:lated);
% regional_MJJ_pdsi_monmean_lon = pdsi_lon_monmean(lonst:loned);
% regional_MJJ_pdsi_monmean_latmat = repmat(...
%     regional_MJJ_pdsi_monmean_lat,[1 length(...
%     regional_MJJ_pdsi_monmean_lon)])';
% regional_MJJ_pdsi_monmean_weighted_mean = nan(1,length(1948:2010))';
% for i = 1:length(1948:2010)
%     regional_MJJ_pdsi_monmean_weighted_mean(i,1) = wmean(...
%         squeeze(regional_MJJ_pdsi_monmean_gridded_mean(:,:,i)),...
%         cosd(regional_MJJ_pdsi_monmean_latmat),'all','omitnan');
% end
% 
% regional_MJJ_pdsi_monmean_weighted_mean_detrended = detrend(regional_MJJ_pdsi_monmean_weighted_mean);
% % Spatially-weighted correlations
% [rins_pdsi_monmean, pins_pdsi_monmean] = corr(detrend(BHIz_ins),regional_MJJ_pdsi_monmean_weighted_mean_detrended,...
%     'Type','Pearson','Rows','pairwise','Tail','right');
% [rrecon_pdsi_monmean, precon_pdsi_monmean] = corr(detrend(RECON(find(YEAR == 1948,1,'first'):find(YEAR == 2009,1,'first'))),...
%     regional_MJJ_pdsi_monmean_weighted_mean_detrended(1:62),'Type','Pearson','Rows','pairwise','Tail','right');
% % Gridded correlations
% [gx,gy,~] = size(regional_MJJ_pdsi_monmean_gridded_mean);
% rins_pdsi_monmean_gridded = nan([gx gy]);
% pins_pdsi_monmean_gridded = nan([gx gy]);
% rrecon_pdsi_monmean_gridded = nan([gx gy]);
% precon_pdsi_monmean_gridded = nan([gx gy]);
% for i = 1:gx
%     for j = 1:gy
%         pdsi_monmean_ts = squeeze(regional_MJJ_pdsi_monmean_gridded_mean(i,j,:));
%         [r,p] = corr(detrend(BHIz_ins),detrend(pdsi_monmean_ts),'Type','Pearson','Rows','pairwise','Tail','right');
%         rins_pdsi_monmean_gridded(i,j) = r;
%         pins_pdsi_monmean_gridded(i,j) = p;
%         [r_r, p_r] = corr(detrend(RECON(find(YEAR == 1948,1,'first'):find(YEAR == 2009,1,'first'))),...
%             detrend(pdsi_monmean_ts(1:62)),'Type','Pearson','Rows','pairwise','Tail','right');
%         rrecon_pdsi_monmean_gridded(i,j) = r_r;
%         precon_pdsi_monmean_gridded(i,j) = p_r;
%     end
% end
% 
% regional_MJJ_precip_monmean_1x1(regional_MJJ_precip_monmean_1x1 < 0) = NaN;
% regional_MJJ_precip_monmean_1x1_gridded_mean = squeeze(mean(regional_MJJ_precip_monmean_1x1,4,'omitnan'));
% latst = find(precip_monmean_1x1_lat <= 38,1,'first');
% lated = find(precip_monmean_1x1_lat >= 26,1,'last');
% lonst = find((precip_monmean_1x1_lon-180) >= 70,1,'first');
% loned = find((precip_monmean_1x1_lon-180) <= 92,1,'last');
% regional_MJJ_precip_monmean_1x1_lat = precip_monmean_1x1_lat(latst:lated);
% regional_MJJ_precip_monmean_1x1_lon = precip_monmean_1x1_lon(lonst:loned);
% regional_MJJ_precip_monmean_1x1_latmat = repmat(...
%     regional_MJJ_precip_monmean_1x1_lat,[1 length(...
%     regional_MJJ_precip_monmean_1x1_lon)])';
% regional_MJJ_precip_monmean_1x1_weighted_mean = nan(1,length(1948:2010))';
% for i = 1:length(1948:2010)
%     regional_MJJ_precip_monmean_1x1_weighted_mean(i,1) = wmean(...
%         squeeze(regional_MJJ_precip_monmean_1x1_gridded_mean(:,:,i)),...
%         cosd(regional_MJJ_precip_monmean_1x1_latmat),'all','omitnan');
% end
% 
% regional_MJJ_precip_monmean_1x1_weighted_mean_detrended = detrend(regional_MJJ_precip_monmean_1x1_weighted_mean);
% % Spatially-weighted correlations
% [rins_precip_monmean_1x1, pins_precip_monmean_1x1] = corr(detrend(BHIz_ins),regional_MJJ_precip_monmean_1x1_weighted_mean_detrended,...
%     'Type','Pearson','Rows','pairwise','Tail','right');
% [rrecon_precip_monmean_1x1, precon_precip_monmean_1x1] = corr(detrend(RECON(find(YEAR == 1948,1,'first'):find(YEAR == 2009,1,'first'))),...
%     regional_MJJ_precip_monmean_1x1_weighted_mean_detrended(1:62),'Type','Pearson','Rows','pairwise','Tail','right');
% % Gridded correlations
% [gx,gy,~] = size(regional_MJJ_precip_monmean_1x1_gridded_mean);
% rins_precip_monmean_1x1_gridded = nan([gx gy]);
% pins_precip_monmean_1x1_gridded = nan([gx gy]);
% rrecon_precip_monmean_1x1_gridded = nan([gx gy]);
% precon_precip_monmean_1x1_gridded = nan([gx gy]);
% for i = 1:gx
%     for j = 1:gy
%         precip_monmean_1x1_ts = squeeze(regional_MJJ_precip_monmean_1x1_gridded_mean(i,j,:));
%         [r,p] = corr(detrend(BHIz_ins),detrend(precip_monmean_1x1_ts),'Type','Pearson','Rows','pairwise','Tail','right');
%         rins_precip_monmean_1x1_gridded(i,j) = r;
%         pins_precip_monmean_1x1_gridded(i,j) = p;
%         [r_r,p_p] = corr(detrend(RECON(find(YEAR == 1948,1,'first'):find(YEAR == 2009,1,'first'))),...
%             detrend(precip_monmean_1x1_ts(1:62)),'Type','Pearson','Rows','pairwise','Tail','right');
%         rrecon_precip_monmean_1x1_gridded(i,j) = r_r;
%         precon_precip_monmean_1x1_gridded(i,j) = p_r;
%     end
% end
% 
% regional_MJJ_soilw_monmean(regional_MJJ_soilw_monmean < 0) = NaN;
% regional_MJJ_soilw_monmean_gridded_mean = squeeze(mean(regional_MJJ_soilw_monmean,4,'omitnan'));
% latst = find(soilw_monmean_lat <= 38,1,'first');
% lated = find(soilw_monmean_lat >= 26,1,'last');
% lonst = find((soilw_monmean_lon-180) >= 70,1,'first');
% loned = find((soilw_monmean_lon-180) <= 92,1,'last');
% regional_MJJ_soilw_monmean_lat = soilw_monmean_lat(latst:lated);
% regional_MJJ_soilw_monmean_lon = soilw_monmean_lon(lonst:loned);
% regional_MJJ_soilw_monmean_latmat = repmat(...
%     regional_MJJ_soilw_monmean_lat,[1 length(...
%     regional_MJJ_soilw_monmean_lon)])';
% regional_MJJ_soilw_monmean_weighted_mean = nan(1,length(1948:2010))';
% for i = 1:length(1948:2010)
%     regional_MJJ_soilw_monmean_weighted_mean(i,1) = wmean(...
%         squeeze(regional_MJJ_soilw_monmean_gridded_mean(:,:,i)),...
%         cosd(regional_MJJ_soilw_monmean_latmat),'all','omitnan');
% end
% 
% % Detrended might not be necessary
% regional_MJJ_soilw_monmean_weighted_mean_detrended = detrend(regional_MJJ_soilw_monmean_weighted_mean);
% % Spatially-weighted correlations
% [rins_soilwd, pins_soilwd] = corr(detrend(BHIz_ins),regional_MJJ_soilw_monmean_weighted_mean_detrended,...
%     'Type','Pearson','Rows','pairwise','Tail','right'); % Spearman gives a slightly higher r-value than Pearson but idk
% [rrecon_soilwd, precon_soilwd] = corr(detrend(RECON(find(YEAR == 1948,1,'first'):find(YEAR == 2009,1,'first'))),...
%     regional_MJJ_soilw_monmean_weighted_mean_detrended(1:62),'Type','Pearson','Rows','pairwise','Tail','right'); % Same thing...
% 
% [rins_soilw, pins_soilw] = corr(BHIz_ins, regional_MJJ_soilw_monmean_weighted_mean,...
%     'Type','Spearman','Rows','complete','Tail','both');
% [rrecon_soilw, precon_soilw] = corr(RECON(find(YEAR == 1948,1,'first'):find(YEAR == 2009,1,'first')),...
%     regional_MJJ_soilw_monmean_weighted_mean,'Type','Spearman','Rows','complete','Tail','both');
% % Gridded correlations
% [gx, gy, ~] = size(regional_MJJ_soilw_monmean_gridded_mean);
% rins_soilw_gridded = nan([gx gy]);
% pins_soilw_gridded = nan([gx gy]);
% rrecon_soilw_gridded = nan([gx gy]);
% precon_soilw_gridded = nan([gx gy]);
% for i = 1:gx
%     for j = 1:gy
%         soilw_ts = squeeze(regional_MJJ_soilw_monmean_gridded_mean(i,j,:));
%         [r,p] = corr(BHIz_ins,soilw_ts,'Type','Spearman','Rows','complete','Tail','both');%detrend(BHIz_ins),detrend(soilw_ts),'Type','Pearson','Rows','pairwise','Tail','right');
%         rins_soilw_gridded(i,j) = r;
%         pins_soilw_gridded(i,j) = p;
%         [r_r, p_r] = corr(RECON(find(YEAR == 1948,1,'first'):find(YEAR == 2009,1,'first')),...
%             soilw_ts(1:62),'Type','Spearman','Rows','complete','Tail','both');
%         %detrend(RECON(find(YEAR == 1948,1,'first'):find(YEAR == 2009,1,'first'))),...
%             %detrend(soilw_ts(1:62)),'Type','Pearson','Rows','pairwise','Tail','right');
%         rrecon_soilw_gridded(i,j) = r_r;
%         precon_soilw_gridded(i,j) = p_r;
%     end
% end
% 
% regional_MJJ_precip_montot_1x1(regional_MJJ_precip_montot_1x1 < 0) = NaN;
% regional_MJJ_precip_montot_1x1_gridded_mean = squeeze(mean(regional_MJJ_precip_montot_1x1,4,'omitnan'));
% latst = find(precip_montot_1x1_lat <= 38,1,'first');
% lated = find(precip_montot_1x1_lat >= 26,1,'last');
% lonst = find((precip_montot_1x1_lon-180) >= 70,1,'first');
% loned = find((precip_montot_1x1_lon-180) <= 92,1,'last');
% regional_MJJ_precip_montot_1x1_lat = precip_montot_1x1_lat(latst:lated);
% regional_MJJ_precip_montot_1x1_lon = precip_montot_1x1_lon(lonst:loned);
% regional_MJJ_precip_montot_1x1_latmat = repmat(...
%     regional_MJJ_precip_montot_1x1_lat,[1 length(...
%     regional_MJJ_precip_montot_1x1_lon)])';
% regional_MJJ_precip_montot_1x1_weighted_mean = nan(1,length(1948:2010))';
% for i = 1:length(1948:2010)
%     regional_MJJ_precip_montot_1x1_weighted_mean(i,1) = wmean(...
%         squeeze(regional_MJJ_precip_montot_1x1_gridded_mean(:,:,i)),...
%         cosd(regional_MJJ_precip_montot_1x1_latmat),'all','omitnan');
% end
% 
% regional_MJJ_precip_montot_1x1_weighted_mean_detrended = detrend(regional_MJJ_precip_montot_1x1_weighted_mean);
% % Spatially-weighted correlations
% [rins_precip_montot_1x1, pins_precip_montot_1x1] = corr(detrend(BHIz_ins),regional_MJJ_precip_montot_1x1_weighted_mean_detrended,...
%     'Type','Pearson','Rows','pairwise','Tail','right');
% [rrecon_precip_montot_1x1, precon_precip_montot_1x1] = corr(detrend(RECON(find(YEAR == 1948,1,'first'):find(YEAR == 2009,1,'first'))),...
%     regional_MJJ_precip_montot_1x1_weighted_mean_detrended(1:62),'Type','Pearson','Rows','pairwise','Tail','right');
% % Gridded correlations
% [gx, gy, ~] = size(regional_MJJ_precip_montot_1x1_gridded_mean);
% rins_precip_montot_1x1_gridded = nan([gx gy]);
% pins_precip_montot_1x1_gridded = nan([gx gy]);
% rrecon_precip_montot_1x1_gridded = nan([gx gy]);
% precon_precip_montot_1x1_gridded = nan([gx gy]);
% for i = 1:gx
%     for j = 1:gy
%         precip_montot_1x1_ts = squeeze(regional_MJJ_precip_montot_1x1_gridded_mean(i,j,:));
%         [r,p] = corr(detrend(BHIz_ins),detrend(precip_montot_1x1_ts),'Type','Pearson','Rows','pairwise','Tail','right');
%         rins_precip_montot_1x1_gridded(i,j) = r;
%         pins_precip_montot_1x1_gridded(i,j) = p;
%         [r_r,p_r] = corr(detrend(RECON(find(YEAR == 1948,1,'first'):find(YEAR == 2009,1,'first'))),...
%             detrend(precip_montot_1x1_ts(1:62)),'Type','Pearson','Rows','pairwise','Tail','right');
%         rrecon_precip_montot_1x1_gridded(i,j) = r_r;
%         precon_precip_montot_1x1_gridded(i,j) = p_r;
%     end
% end
% 

%
%% Correlations with instrumental records

% July 14, 2022

% [rmjjpdsicruwmean_, rmjjpdsicruwmean_p] = corr(BHIz_ins,...
%     regional_MJJ_pdsi_cru_weighted_mean,'Type','Spearman',...
%     'Tail','both','Rows','complete');
% [rmjjpdsimmeanwmean_r, rmjjpdsimmeanwmean_p] = corr(BHIz_ins,...
%     regional_MJJ_pdsi_monmean_weighted_mean,'Type','Spearman',...
%     'Tail','both','Rows','complete');
% 
% 
% figure
% 
% scatter(BHIz_ins,regional_MJJ_pdsi_cru_weighted_mean)
% xlabel('BHI')
% ylabel('Mean MJJ PDSI (CRU)')
% 
% scatter(BHIz_ins, regional_MJJ_pdsi_monmean_weighted_mean)
% xlabel('BHI')
% ylabel('Monthly Mean MJJ PDSI')
% 
% [plon_m_soilw,plat_m_soilw] = meshgrid(regional_MJJ_soilw_monmean_lon-360,regional_MJJ_soilw_monmean_lat);
% figure
% subplot(2,2,[1 2])
% pcolor(plon_m_soilw,plat_m_soilw,rins_soilw_gridded')
% shading interp
% usbasemap
% cmocean('diff')
% xlim([-107 -70])
% ylim([23 40])
% set(gca,'Layer','top')
% set(gca,'clim',[-1 1]);
% hc = colorbar('vert');
% set(hc, 'YAxisLocation','right')
% hc.Label.String = 'Correlation coefficient';
% hc.Label.FontSize = 18;
% hc.Label.FontName = 'Helvetica';
% hc.Ticks = (-1:0.1:1);
% hc.TickLabels = (-1:0.2:1);
% hold on
% maskp = pins_soilw_gridded < 0.05;
% stipple(plon_m_soilw,plat_m_soilw,maskp','density',250,'marker','x');
% v = [0.05,0.05];
% contour(plon_m_soilw,plat_m_soilw,pins_soilw_gridded',v,'LineColor','k');
% hold off

%
%% BHI Recon and Other Recons

% % Load various datasets
% % Luterbacher et al. 2000 NAO recon
%     % Removing the non-summer values
% % LseaNAOind = Luterbacheretal2001seaNAO15001658.NAO; 
% % LseaNAOyr = Luterbacheretal2001seaNAO15001658.Year;
% % LseaNAOseanum = Luterbacheretal2001seaNAO15001658.SEAnum;
% % LseaNAOind(LseaNAOseanum < 4) = [];
% % LseaNAOyr(LseaNAOseanum < 4) = [];
% % 
% % LmonNAOind = Luterbacheretal2001monNAO.NAO;
% % LmonNAOyr = Luterbacheretal2001monNAO.Year;
% % LmonNAOmonnum = Luterbacheretal2001monNAO.MonNum;
% % LmonNAOind(LmonNAOmonnum < 6) = [];
% % LmonNAOyr(LmonNAOmonnum < 6) = [];
% % LmonNAOmonnum(LmonNAOmonnum < 6) = [];
% % LmonNAOind(LmonNAOmonnum > 8) = [];
% % LmonNAOyr(LmonNAOmonnum > 8) = [];
% % LmonNAOmonnum(LmonNAOmonnum > 8) = [];
% % LmonNAOind = reshape(LmonNAOind,[3,342]);
% % LmonNAOindJJAave = mean(LmonNAOind,1);
% % LmonNAOindJJAave = LmonNAOindJJAave';
% LmonNAOyrJJAave = unique(LmonNAOyr);
% LNAOyr = [LseaNAOyr; LmonNAOyrJJAave];
% LNAOind = [LseaNAOind; LmonNAOindJJAave];

%
%% Plotting BHI Against Other Recons

% % BHI: 1430–2011
% % Hadley Position: 1430–2003
% % Summer NAO (Folland): 1706–2009
% % TSI: 1610–2000
% % Summer NAO (Luterbacher): 1500–2000
% % WAtl SST: 1552–2007
% % US SAT: 1430–1987
% 
% % For now, let's cut it off at 2000. 
% % BHI -> TSI -> WAtl SST -> US SAT -> Summer NAO (both) -> Hadley
% % Add mean line later
% BHIedy_loc = find(YEAR == 2000,1,'first');
% 
% LeanTSIyr = floor(Lean2000TSI.YEAR1);
% LeanTSIyr = [(1430:1:1609)'; LeanTSIyr];
% TSInans = nan(length(LeanTSIyr(1:180)));
% LeanTSI11cy = Lean2000TSI.yrCYCLE11yrC;
% LeanTSI11cy = [TSInans(1:180,1); LeanTSI11cy];
% LeanTSIbkgd = Lean2000TSI.YCLEBKGRND;
% LeanTSIbkgd = [TSInans(1:180,1); LeanTSIbkgd];
% deeporange = rgb('deep orange');
% black = rgb('black');
% 
% WAtl_yr = Tierneyetal2015WAtl.age_AD;
% WAtl_yr = flipud(WAtl_yr);
% WAtl_yr_loc = find(WAtl_yr == 2000,1,'first');
% WAtl_yr = WAtl_yr(1:WAtl_yr_loc);
% WAtl_SSTA = Tierneyetal2015WAtl.sstaprmaranom;
% WAtl_SSTA = flipud(WAtl_SSTA(1:WAtl_yr_loc));
% WAtl_yr = [(1430:1:1551)'; WAtl_yr];
% WAtl_nans = nan(length(1430:1:1551),1);
% WAtl_SSTA = [WAtl_nans; WAtl_SSTA];
% coral = rgb('coral');
% 
% SATUS_yr = Trouetetal2013SATUS.age_AD;
% SATUS_temp = Trouetetal2013SATUS.temp2mann;
% SATUS_yr = [SATUS_yr; (1988:1:2000)'];
% SATUS_nans = nan(length(1988:1:2000),1);
% SATUS_temp = [SATUS_temp; SATUS_nans];
% brownishred = rgb('brownish red');
% 
% LNAOyr = [LseaNAOyr;LmonNAOyr];
% LNAOind = [LseaNAOind; LmonNAOindJJAave];
% LNAOyr = [(1430:1:1499)';LNAOyr];
% LNAO_nans = nan(length(1430:1:1499),1);
% LNAOind = [LNAO_nans; LNAOind];
% kellygreen = rgb('kelly green');
% FNAOyr = Follandetal2009summerNAO.Year;
% FNAOyr_loc = find(FNAOyr == 2000,1,'first');
% FNAOyr = FNAOyr(1:FNAOyr_loc);
% FNAOind = Follandetal2009summerNAO.SNAOrec;
% FNAOind = FNAOind(1:FNAOyr_loc);
% FNAOyr = [(1430:1:1705)'; FNAOyr];
% FNAO_nans = nan(length(1430:1:1705),1);
% FNAOind = [FNAO_nans; FNAOind];
% ochre = rgb('ochre');
% 
% Hadyr = AlfaroSanchezetal2018HadleyPosition.age_CE;
% Hadlat = AlfaroSanchezetal2018HadleyPosition.TBNB;
% Hadyr_loc = find(Hadyr == 2000,1,'first');
% Hadyr = Hadyr(1:Hadyr_loc);
% Hadlat = Hadlat(1:Hadyr_loc);
% deeplavender = rgb('deep lavender');
% 
% panel = {'a','b','c','d','e','f'};
% figure
% ax(1) = subsubplot(6,1,1);
% plot(YEAR(1:BHIedy_loc),RECON(1:BHIedy_loc),'-','Color',waterblue,'LineWidth',2)
% hold on
% patch([YEAR(1:BHIedy_loc); flipud(YEAR(1:BHIedy_loc))],[PREUP_95(1:BHIedy_loc); flipud(PRELO_05(1:BHIedy_loc))],waterblue,'FaceAlpha',0.3,'LineStyle','none')
% hold off
% axis tight
% ntitle(panel{1},'location','nw')
% ax(2) = subsubplot(6,1,2);
% %plot(LeanTSIyr,LeanTSI11cy,'-','Color',deeporange,'LineWidth',2)
% %hold on
% plot(YEAR(1:BHIedy_loc),LeanTSIbkgd,'-','Color',deeporange,'LineWidth',2)
% %hold off
% axis tight
% ntitle(panel{2},'location','nw')
% ax(3) = subsubplot(6,1,3);
% plot(YEAR(1:BHIedy_loc),WAtl_SSTA,'-','Color',coral,'LineWidth',2)
% axis tight
% ntitle(panel{3},'location','nw')
% ax(4) = subsubplot(6,1,4);
% plot(YEAR(1:BHIedy_loc),SATUS_temp,'-','Color',brownishred,'LineWidth',2);
% axis tight
% ntitle(panel{4},'location','nw')
% ax(5) = subsubplot(6,1,5);
% plot(YEAR(1:BHIedy_loc), LNAOind,'-','Color',kellygreen,'LineWidth',2);
% hold on
% plot(YEAR(1:BHIedy_loc),FNAOind,'-','Color',ochre,'LineWidth',2);
% hold off
% axis tight
% ntitle(panel{5},'location','nw')
% ax(6) = subsubplot(6,1,6);
% plot(YEAR(1:BHIedy_loc),Hadlat,'-','Color',deeplavender,'LineWidth',2)
% axis tight
% ntitle(panel{6},'location','nw')
% xlabel('Year (CE)')
% 
% set(ax(2:2:6),'YAxisLocation','right')
% 

%
%% Try this for the bit above

% panel = {'A','B','C','D','E','F','G'};
% SATUS25bound = Trouetetal2013SATUS.temp2mann25;
% SATUS975bound = Trouetetal2013SATUS.temp2mann975;
% Hadleylat_se = AlfaroSanchezetal2018HadleyPosition.TBNB1se;
% Hadleylat_se1 = AlfaroSanchezetal2018HadleyPosition.TBNB1se1;
% %figure
% t = tiledlayout(7,1,'TileSpacing','none');
% ax1 = nexttile;%subplot(6,2,[1,2])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc),RECON(1:BHIedy_loc),'-','Color',waterblue,'LineWidth',2)
% hold on
% patch([YEAR(1:BHIedy_loc); flipud(YEAR(1:BHIedy_loc))],[PREUP_95(1:BHIedy_loc); flipud(PRELO_05(1:BHIedy_loc))],waterblue,'FaceAlpha',0.3,'LineStyle','none')
% hold off
% axis tight
% ntitle(panel{1},'location','nw')
% ax1.Box = 'off';
% %ax.XColor = 'none';
% ax1.XAxisLocation = 'top';
% ntitle(panel{1},'location','nw')
% ylabel('BHI')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax2 = nexttile;%subplot(6,2,[3,4])
% anomaly(YEAR,RECON,'thresh',mean(BHIz_ins,'omitnan'))
% axis tight
% ax2.Box = 'off';
% ax2.YColor = 'k';
% ax2.XColor = 'none';
% xlim([1430 2000]);
% ax2.YAxisLocation = 'right';
% ylabel('BHI Anomaly')
% ntitle(panel{2},'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax3 = nexttile;
% %ax = gca;
% yyaxis left
% h = bar(age_AD,SO4_flux,'FaceColor',rgb('gunmetal'),'EdgeColor',rgb('gunmetal'))
% %h
% ylabel('SO_2_- flux (kg km^-^2)')
% h.ShowBaseLine = 'off';%ax2.XColor = 'none';
% h.BarWidth = 1;
% ax3.YColor = rgb('gunmetal');%'k';
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% yyaxis right
% plot(YEAR(1:BHIedy_loc),LeanTSIbkgd,'-','Color',deeporange,'LineWidth',2)
% ax3.YColor = deeporange;%'k';
% axis tight
% ax3.Box = 'off';
% ax3.XColor = 'none';
% xlim([1430 2000]);
% %ax2.YAxisLocation = 'right';
% ylabel('TSI (W m^-^2)')
% % yyaxis left
% % bar(age_AD,SO4_flux,'FaceColor',rgb('gunmetal'),'EdgeColor',rgb('gunmetal'))
% % %h
% % ylabel('SO_2_- flux (kg km^-^2)')
% % ax2.XColor = 'none';
% % ax2.YColor = 'k';
% %hold off
% ntitle(panel{3}, 'location','nw')
% %vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax4 = nexttile;%subplot(6,2,[5,6])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc),WAtl_SSTA,'-','Color',coral,'LineWidth',2)
% axis tight
% ax4.Box = 'off';
% ax4.XColor = 'none';
% xlim([1430 2000]);
% ax4.YAxisLocation = 'right';
% ylabel('SSTA (^oC)')
% ntitle(panel{4}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax5 = nexttile;%subplot(6,2,[7,8])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc),SATUS_temp,'-','Color',brownishred,'LineWidth',2);
% hold on
% patch([YEAR(1:558); flipud(YEAR(1:558))],[SATUS975bound; flipud(SATUS25bound)],brownishred,'FaceAlpha',0.3,'LineStyle','none')
% hold off
% axis tight
% ax5.Box = 'off';
% ax5.XColor = 'none';
% xlim([1430 2000])
% %ax5.YAxisLocation = 'right';
% ylabel('SATA (^oC)')
% ntitle(panel{5}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax6 = nexttile;%subplot(6,2,[9,10])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc), LNAOind,'-','Color',kellygreen,'LineWidth',2);
% % hold on
% % plot(YEAR(1:BHIedy_loc),FNAOind,'-','Color',ochre,'LineWidth',2);
% % hold off
% ax6.Box = 'off';
% ax6.XColor = 'none';
% xlim([1430 2000])
% ylabel('Summer NAO')
% ax6.YAxisLocation = 'right';
% ntitle(panel{6}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax7 = nexttile;%subplot(6,2,[11,12])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc),Hadlat,'-','Color',deeplavender,'LineWidth',2)
% hold on
% patch([YEAR(1:BHIedy_loc); flipud(YEAR(1:BHIedy_loc))],[Hadleylat_se1(1:BHIedy_loc); flipud(Hadleylat_se(1:BHIedy_loc))],deeplavender,'FaceAlpha',0.3,'LineStyle','none')
% hold off
% axis tight
% ax7.Box = 'off';
% xlim([1430 2000])
% %ax7.YAxisLocation = 'right';
% ylabel('Hadley Latitude (^oN)')
% xlabel('Years (CE)')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ntitle(panel{7},'location','nw')
% 

%
%% Same for the above section but with 5-yr moving averages

% WAtl_SSTA_5mm = movmean(WAtl_SSTA,5,'omitnan','endpoints','fill');
% SATUS_temp_5mm = movmean(SATUS_temp,5,'omitnan','endpoints','fill');
% SATUS_25bound_5mm = movmean(SATUS25bound,5,'omitnan','endpoints','fill');
% SATUS_975bound_5mm = movmean(SATUS975bound,5,'omitnan','endpoints','fill');
% LNAOind_5mm = movmean(LNAOind,5,'omitnan','endpoints','fill');
% Hadlat_5mm = movmean(Hadlat,5,'omitnan','endpoints','fill');
% Hadlat_se_5mm = movmean(Hadleylat_se,5,'omitnan','endpoints','fill');
% Hadlat_se1_5mm = movmean(Hadleylat_se1,5,'omitnan','endpoints','fill');

%
%% Plotting the section above

% panel = {'A','B','C','D','E','F','G'};
% t2 = tiledlayout(7,1,'TileSpacing','none');
% ax12 = nexttile;%subplot(6,2,[1,2])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc),BHIz_recon_5movmean(1:BHIedy_loc),'-','Color',waterblue,'LineWidth',2)
% hold on
% %patch([YEAR(1:BHIedy_loc); flipud(YEAR(1:BHIedy_loc))],[BHIz_recon_pup95_5movmean(1:BHIedy_loc); flipud(BHIz_recon_plo05_5movmean(1:BHIedy_loc))],waterblue,'FaceAlpha',0.3,'LineStyle','none')
% patch([YEAR(3:580); flipud(YEAR(3:580))],[BHIz_recon_pup95_5movmean(3:580); flipud(BHIz_recon_plo05_5movmean(3:580))],waterblue,'FaceAlpha',0.3,'LineStyle','none')
% hold off
% axis tight
% ntitle(panel{1},'location','nw')
% ax12.Box = 'off';
% %ax.XColor = 'none';
% ax12.XAxisLocation = 'top';
% ntitle(panel{1},'location','nw')
% ylabel('BHIz')
% xlim([1430 2000])
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax22 = nexttile;%subplot(6,2,[3,4])
% anomaly(YEAR,BHIz_recon_5movmean,'thresh',mean(BHIz_ins_5mm,'omitnan'))%anomaly(YEAR,RECON,'thresh',mean(BHIz_ins,'omitnan'))
% axis tight
% ax22.Box = 'off';
% ax22.YColor = 'k';
% ax22.XColor = 'none';
% xlim([1430 2000]);
% ax22.YAxisLocation = 'right';
% ylabel('BHI Anomaly (5-yr mean)')
% ntitle(panel{2},'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax32 = nexttile;
% yyaxis left
% h32 = bar(age_AD,SO4_flux,'FaceColor',rgb('gunmetal'),'EdgeColor',rgb('gunmetal'))
% %h
% ylabel('SO_2_- flux (kg km^-^2)')
% h32.ShowBaseLine = 'off';%ax2.XColor = 'none';
% h32.BarWidth = 3;
% ax32.YColor = rgb('gunmetal');%'k';
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% %ax = gca;
% yyaxis right
% plot(YEAR(1:BHIedy_loc),LeanTSIbkgd,'-','Color',deeporange,'LineWidth',2)
% axis tight
% ax32.Box = 'off';
% ax32.XColor = 'none';
% ax32.YColor = deeporange;%'k';
% xlim([1430 2000]);
% %ax22.YAxisLocation = 'right';
% ylabel('TSI (W m^-^2)')
% ntitle(panel{3}, 'location','nw')
% %vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax42 = nexttile;%subplot(6,2,[5,6])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc),WAtl_SSTA_5mm,'-','Color',coral,'LineWidth',2)
% axis tight
% ax42.Box = 'off';
% ax42.XColor = 'none';
% ax42.YAxisLocation = 'right';
% xlim([1430 2000]);
% ylabel('SSTA (^oC; 5-yr mean)')
% ntitle(panel{4}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax52 = nexttile;%subplot(6,2,[7,8])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc),SATUS_temp_5mm,'-','Color',brownishred,'LineWidth',2);
% hold on
% patch([YEAR(1:558); flipud(YEAR(1:558))],[SATUS_975bound_5mm; flipud(SATUS_25bound_5mm)],brownishred,'FaceAlpha',0.3,'LineStyle','none')
% hold off
% axis tight
% ax52.Box = 'off';
% ax52.XColor = 'none';
% xlim([1430 2000])
% %ax52.YAxisLocation = 'right';
% ylabel('SATA (^oC; 5-yr mean)')
% ntitle(panel{5}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax62 = nexttile;%subplot(6,2,[9,10])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc), LNAOind_5mm,'-','Color',kellygreen,'LineWidth',2);
% % hold on
% % plot(YEAR(1:BHIedy_loc),FNAOind,'-','Color',ochre,'LineWidth',2);
% % hold off
% ax62.Box = 'off';
% ax62.XColor = 'none';
% ax62.YAxisLocation = 'right';
% xlim([1430 2000])
% ylabel('Summer NAO (5-yr mean)')
% ntitle(panel{6}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax72 = nexttile;%subplot(6,2,[11,12])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc),Hadlat_5mm,'-','Color',deeplavender,'LineWidth',2)
% hold on
% patch([YEAR(1:BHIedy_loc); flipud(YEAR(1:BHIedy_loc))],[Hadlat_se1_5mm(1:BHIedy_loc); flipud(Hadlat_se_5mm(1:BHIedy_loc))],deeplavender,'FaceAlpha',0.3,'LineStyle','none')
% hold off
% axis tight
% ax72.Box = 'off';
% xlim([1430 2000])
% %ax72.YAxisLocation = 'right';
% ylabel('Hadley Latitude (^oN; 5-yr mean)')
% xlabel('Years (CE)')
% ntitle(panel{7},'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% 

%
%% Same thing but with 10-yr moving averages

% WAtl_SSTA_10mm = movmean(WAtl_SSTA,10,'omitnan','endpoints','fill');
% SATUS_temp_10mm = movmean(SATUS_temp,10,'omitnan','endpoints','fill');
% SATUS_25bound_10mm = movmean(SATUS25bound,10,'omitnan','endpoints','fill');
% SATUS_975bound_10mm = movmean(SATUS975bound,10,'omitnan','endpoints','fill');
% LNAOind_10mm = movmean(LNAOind,10,'omitnan','endpoints','fill');
% Hadlat_10mm = movmean(Hadlat,10,'omitnan','endpoints','fill');
% Hadlat_se_10mm = movmean(Hadleylat_se,10,'omitnan','endpoints','fill');
% Hadlat_se1_10mm = movmean(Hadleylat_se1,10,'omitnan','endpoints','fill');
% 
% t3 = tiledlayout(6,1,'TileSpacing','none');
% ax13 = nexttile;%subplot(6,2,[1,2])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc),BHIz_recon_10movmean(1:BHIedy_loc),'-','Color',waterblue,'LineWidth',2)
% hold on
% patch([YEAR(1:BHIedy_loc); flipud(YEAR(1:BHIedy_loc))],[PREUP_95(1:BHIedy_loc); flipud(PRELO_05(1:BHIedy_loc))],waterblue,'FaceAlpha',0.3,'LineStyle','none')
% hold off
% axis tight
% ntitle(panel{1},'location','nw')
% ax12.Box = 'off';
% %ax.XColor = 'none';
% ax12.XAxisLocation = 'top';
% ntitle(panel{1},'location','nw')
% ylabel('BHIz')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax22 = nexttile;%subplot(6,2,[3,4])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc),LeanTSIbkgd,'-','Color',deeporange,'LineWidth',2)
% axis tight
% ax22.Box = 'off';
% ax22.XColor = 'none';
% xlim([1430 2000]);
% ax22.YAxisLocation = 'right';
% ylabel('TSI')
% ntitle(panel{2}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax32 = nexttile;%subplot(6,2,[5,6])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc),WAtl_SSTA_5mm,'-','Color',coral,'LineWidth',2)
% axis tight
% ax32.Box = 'off';
% ax32.XColor = 'none';
% xlim([1430 2000]);
% ylabel('SSTA (^oC)')
% ntitle(panel{3}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax42 = nexttile;%subplot(6,2,[7,8])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc),SATUS_temp_5mm,'-','Color',brownishred,'LineWidth',2);
% hold on
% patch([YEAR(1:558); flipud(YEAR(1:558))],[SATUS_975bound_5mm; flipud(SATUS_25bound_5mm)],brownishred,'FaceAlpha',0.3,'LineStyle','none')
% hold off
% axis tight
% ax42.Box = 'off';
% ax42.XColor = 'none';
% xlim([1430 2000])
% ax42.YAxisLocation = 'right';
% ylabel('SATA (^oC)')
% ntitle(panel{4}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax52 = nexttile;%subplot(6,2,[9,10])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc), LNAOind_5mm,'-','Color',kellygreen,'LineWidth',2);
% % hold on
% % plot(YEAR(1:BHIedy_loc),FNAOind,'-','Color',ochre,'LineWidth',2);
% % hold off
% ax52.Box = 'off';
% ax52.XColor = 'none';
% xlim([1430 2000])
% ylabel('Summer NAO')
% ntitle(panel{5}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax62 = nexttile;%subplot(6,2,[11,12])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc),Hadlat_5mm,'-','Color',deeplavender,'LineWidth',2)
% hold on
% patch([YEAR(1:BHIedy_loc); flipud(YEAR(1:BHIedy_loc))],[Hadlat_se1_5mm(1:BHIedy_loc); flipud(Hadlat_se_5mm(1:BHIedy_loc))],deeplavender,'FaceAlpha',0.3,'LineStyle','none')
% hold off
% axis tight
% ax62.Box = 'off';
% xlim([1430 2000])
% ax62.YAxisLocation = 'right';
% ylabel('Hadley Latitude (^oN)')
% xlabel('Years (CE)')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');

%
%% Annual Anomaly Plots

% panel = {'1','2','3','4','5','6'};
% t4 = tiledlayout(6,1,'TileSpacing','none');
% ax14 = nexttile;%subplot(6,2,[1,2])
% %ax = gca;
% %plot(YEAR(1:BHIedy_loc),RECON(1:BHIedy_loc),'-','Color',waterblue,'LineWidth',2)
% %anomaly(YEAR,RECON)%YEAR(1:BHIedy_loc),RECON(1:BHIedy_loc))
% anomaly(YEAR,BHIz_recon_5movmean,'thresh',mean(BHIz_recon_5movmean,'omitnan'))
% %hold on
% %patch([YEAR(1:BHIedy_loc); flipud(YEAR(1:BHIedy_loc))],[PREUP_95(1:BHIedy_loc); flipud(PRELO_05(1:BHIedy_loc))],waterblue,'FaceAlpha',0.3,'LineStyle','none')
% %hold off
% axis tight
% ntitle(panel{1},'location','nw')
% xlim([1430 2000])
% ax14.Box = 'off';
% %ax.XColor = 'none';
% ax14.XAxisLocation = 'top';
% ax14.XColor = 'k';
% ntitle(panel{1},'location','nw')
% ylabel('BHI Anomaly (5-yr mean)')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax24 = nexttile;%subplot(6,2,[3,4])
% %ax = gca;
% plot(YEAR(1:BHIedy_loc),LeanTSIbkgd,'-','Color',deeporange,'LineWidth',2)
% axis tight
% ax24.Box = 'off';
% ax24.XColor = 'none';
% xlim([1430 2000]);
% ax24.YAxisLocation = 'right';
% ylabel('TSI (W m^-^2; 5-yr mean)')
% ntitle(panel{2}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax34 = nexttile;%subplot(6,2,[5,6])
% %ax = gca;
% %plot(YEAR(1:BHIedy_loc),WAtl_SSTA,'-','Color',coral,'LineWidth',2)
% plot(YEAR(1:BHIedy_loc),WAtl_SSTA_5mm,'-','Color',coral,'LineWidth',2)
% axis tight
% ax34.Box = 'off';
% ax34.XColor = 'none';
% xlim([1430 2000]);
% ylabel('SSTA (^oC; 5-yr mean)')
% ntitle(panel{3}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax44 = nexttile;%subplot(6,2,[7,8])
% %ax = gca;
% %plot(YEAR(1:BHIedy_loc),SATUS_temp,'-','Color',brownishred,'LineWidth',2);
% plot(YEAR(1:BHIedy_loc),SATUS_temp_5mm,'-','Color',brownishred,'LineWidth',2);
% hold on
% patch([YEAR(1:558); flipud(YEAR(1:558))],[SATUS_975bound_5mm; flipud(SATUS_25bound_5mm)],brownishred,'FaceAlpha',0.3,'LineStyle','none')
% %patch([YEAR(1:558); flipud(YEAR(1:558))],[SATUS975bound; flipud(SATUS25bound)],brownishred,'FaceAlpha',0.3,'LineStyle','none')
% hold off
% axis tight
% ax44.Box = 'off';
% ax44.XColor = 'none';
% xlim([1430 2000])
% ax44.YAxisLocation = 'right';
% ylabel('SATA (^oC; 5-yr mean)')
% ntitle(panel{4}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax54 = nexttile;%subplot(6,2,[9,10])
% %ax = gca;
% %anomaly(YEAR(1:BHIedy_loc), LNAOind_5mm);%YEAR(1:BHIedy_loc), LNAOind);
% plot(YEAR(1:BHIedy_loc), LNAOind_5mm,'-','Color',kellygreen,'LineWidth',2)%LNAOind,'-','Color',kellygreen,'LineWidth',2);
% % hold on
% % plot(YEAR(1:BHIedy_loc),FNAOind,'-','Color',ochre,'LineWidth',2);
% % hold off
% ax54.Box = 'off';
% ax54.XColor = 'none';
% xlim([1430 2000])
% ylabel('Summer NAO Anomaly (5-yr mean)')
% ntitle(panel{5}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax64 = nexttile;%subplot(6,2,[11,12])
% %ax = gca;
% %plot(YEAR(1:BHIedy_loc),Hadlat,'-','Color',deeplavender,'LineWidth',2)
% anomaly(YEAR(1:BHIedy_loc),Hadlat_5mm,'thresh',mean(Hadlat_5mm,'omitnan'));%YEAR(1:BHIedy_loc),Hadlat,'thresh',mean(Hadlat,'omitnan'))
% %hold on
% %patch([YEAR(1:BHIedy_loc); flipud(YEAR(1:BHIedy_loc))],[Hadleylat_se1(1:BHIedy_loc); flipud(Hadleylat_se(1:BHIedy_loc))],deeplavender,'FaceAlpha',0.3,'LineStyle','none')
% %hold off
% axis tight
% ax64.Box = 'off';
% xlim([1430 2000])
% ax64.YAxisLocation = 'right';
% ylabel('Hadley Latitude Anomaly (^oN; 5-yr mean)')
% xlabel('Years (CE)')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ntitle(panel{6},'location','nw')

%
%% Guess we're plotting some Gaussian filters

% % Nixing the anomaly plot and the Summer NAO plot
% bhirecon_20ygauss = smoothdata(RECON,'gaussian',20);
% bhirecon_preup95_20ygauss = smoothdata(PREUP_95,'gaussian',20);
% bhirecon_prelo05_20ygauss = smoothdata(PRELO_05,'gaussian',20);
% watl_ssta_20ygauss = smoothdata(WAtl_SSTA,'gaussian',20);
% satus_temp_20ygauss = smoothdata(SATUS_temp(1:558),'gaussian',20);
% satus_975bound_20ygauss = smoothdata(SATUS975bound,'gaussian',20);
% satus_25bound_20ygauss = smoothdata(SATUS25bound,'gaussian',20);
% hadlat_20ygauss = smoothdata(Hadlat,'gaussian',20);
% hadlat_se_20ygauss = smoothdata(Hadlat_se(1:571),'gaussian',20);
% hadlat_se1_20ygauss = smoothdata(Hadlat_se1(1:571),'gaussian',20);
% GeoB9501_dustprop_20ygauss = smoothdata(GeoB9501_dustprop,'gaussian',20);
% GeoB9501_dustmin_20ygauss = smoothdata(GeoB9501_dustmin,'gaussian',20);
% GeoB9501_dustmax_20ygauss = smoothdata(GeoB9501_dustmax,'gaussian',20);

%
%% Okay now I'm plotting


% panel = {'A','B','C','D','E','F','G'};
% t5 = tiledlayout(7,1,'TileSpacing','none');
% ax15 = nexttile;%subplot(6,2,[1,2])
% %ax = gca;
% %plot(YEAR(1:BHIedy_loc),BHIz_recon_5movmean(1:BHIedy_loc),'-','Color',waterblue,'LineWidth',2)
% plot(YEAR,bhirecon_20ygauss,'-','Color',waterblue,'LineWidth',3)
% hold on
% %patch([YEAR(1:BHIedy_loc); flipud(YEAR(1:BHIedy_loc))],[BHIz_recon_pup95_5movmean(1:BHIedy_loc); flipud(BHIz_recon_plo05_5movmean(1:BHIedy_loc))],waterblue,'FaceAlpha',0.3,'LineStyle','none')
% %patch([YEAR(3:580); flipud(YEAR(3:580))],[BHIz_recon_pup95_5movmean(3:580); flipud(BHIz_recon_plo05_5movmean(3:580))],waterblue,'FaceAlpha',0.3,'LineStyle','none')
% patch([YEAR; flipud(YEAR)],[bhirecon_preup95_20ygauss; flipud(bhirecon_prelo05_20ygauss)],waterblue,'FaceAlpha',0.3,'LineStyle','none')
% hold off
% axis tight
% ntitle(panel{1},'location','nw')
% ax15.Box = 'off';
% %ax.XColor = 'none';
% ax15.XAxisLocation = 'top';
% ax15.YAxisLocation = 'right';
% ntitle(panel{1},'location','nw')
% ylabel('BHI')
% xlim([1430 2000])
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% vfill(1797,1820,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax25 = nexttile;%subplot(6,2,[3,4])
% %anomaly(YEAR,BHIz_recon_5movmean,'thresh',mean(BHIz_ins_5mm,'omitnan'))%anomaly(YEAR,RECON,'thresh',mean(BHIz_ins,'omitnan'))
% % axis tight
% % ax22.Box = 'off';
% % ax22.YColor = 'k';
% % ax22.XColor = 'none';
% % xlim([1430 2000]);
% % ax22.YAxisLocation = 'right';
% % ylabel('BHI Anomaly (5-yr mean)')
% % ntitle(panel{2},'location','nw')
% % vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% % ax32 = nexttile;
% yyaxis left
% h25 = bar(age_AD,SO4_flux,'FaceColor',rgb('gunmetal'),'EdgeColor',rgb('gunmetal'))
% %h
% ylabel('SO_2_- flux (kg km^-^2)')
% h25.ShowBaseLine = 'off';%ax2.XColor = 'none';
% h25.BarWidth = 3;
% ax25.YColor = rgb('gunmetal');%'k';
% %ax25.YTick = 0:10:40;
% ax25.YMinorTick = 'on';
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% vfill(1797,1820,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% %ax = gca;
% yyaxis right
% plot(YEAR(1:BHIedy_loc),LeanTSIbkgd,'-','Color',deeporange,'LineWidth',2)
% axis tight
% ax25.Box = 'off';
% ax25.XColor = 'none';
% ax25.YColor = deeporange;%'k';
% xlim([1430 2000]);
% %ax22.YAxisLocation = 'right';
% ylabel('TSI (W m^-^2)')
% ntitle(panel{2}, 'location','nw')
% %vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax35 = nexttile;%subplot(6,2,[5,6])
% %ax = gca;
% %plot(YEAR(1:BHIedy_loc),WAtl_SSTA_5mm,'-','Color',coral,'LineWidth',2)
% plot(WAtl_yr, watl_ssta_20ygauss,'-','Color',coral,'LineWidth',3)
% axis tight
% ax35.Box = 'off';
% ax35.XColor = 'none';
% ax35.YAxisLocation = 'right';
% xlim([1430 2000]);
% ylabel('SSTA (^oC)')
% ntitle(panel{3}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% vfill(1797,1820,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax45 = nexttile;%subplot(6,2,[7,8])
% %ax = gca;
% %plot(YEAR(1:BHIedy_loc),SATUS_temp_5mm,'-','Color',brownishred,'LineWidth',2);
% plot(SATUS_yr(1:558),satus_temp_20ygauss,'-','Color',brownishred,'LineWidth',3);
% hold on
% %patch([YEAR(1:558); flipud(YEAR(1:558))],[SATUS_975bound_5mm; flipud(SATUS_25bound_5mm)],brownishred,'FaceAlpha',0.3,'LineStyle','none')
% patch([SATUS_yr(1:558); flipud(SATUS_yr(1:558))],[satus_975bound_20ygauss; flipud(satus_25bound_20ygauss)],brownishred,'FaceAlpha',0.3,'LineStyle','none')
% hold off
% % hold on
% % indigo = rgb('indigo');
% % plot(SCSAT_year,SCSAT_20ygauss,'-','Color',indigo,'LineWidth',3)
% % hold off
% % hold on
% % sunyellow = rgb('sun yellow');
% % plot(FLSAT_year,FLSAT_20ygauss,'-','Color',sunyellow,'LineWidth',3)
% % hold off
% axis tight
% ax45.Box = 'off';
% ax45.XColor = 'none';
% %ax45.YAxisLocation = 'right';
% xlim([1430 2000])
% %ax52.YAxisLocation = 'right';
% ylabel('SATA (^oC)')
% ntitle(panel{4}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% vfill(1797,1820,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% % ax62 = nexttile;%subplot(6,2,[9,10])
% % %ax = gca;
% % plot(YEAR(1:BHIedy_loc), LNAOind_5mm,'-','Color',kellygreen,'LineWidth',2);
% % % hold on
% % % plot(YEAR(1:BHIedy_loc),FNAOind,'-','Color',ochre,'LineWidth',2);
% % % hold off
% % ax62.Box = 'off';
% % ax62.XColor = 'none';
% % ax62.YAxisLocation = 'right';
% % xlim([1430 2000])
% % ylabel('Summer NAO (5-yr mean)')
% % ntitle(panel{6}, 'location','nw')
% % vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax55 = nexttile;%subplot(6,2,[11,12])
% indigo = rgb('indigo');
% plot(SCSAT_year,SCSAT_20ygauss,'-','Color',indigo,'LineWidth',3)
% %hold off
% hold on
% sunyellow = rgb('sun yellow');
% plot(FLSAT_year,FLSAT_20ygauss,'-','Color',sunyellow,'LineWidth',3)
% hold off
% axis tight
% ax55.Box = 'off';
% ax55.XColor = 'none';
% ax55.YAxisLocation = 'right';
% xlim([1430 2000])
% ylabel('SATA (^oC)')
% ntitle(panel{5}, 'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% vfill(1797,1820,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% %ax = gca;
% %plot(YEAR(1:BHIedy_loc),Hadlat_5mm,'-','Color',deeplavender,'LineWidth',2)
% ax65 = nexttile;
% plot(Hadyr,hadlat_20ygauss,'-','Color',deeplavender,'LineWidth',3)
% hold on
% %patch([YEAR(1:BHIedy_loc); flipud(YEAR(1:BHIedy_loc))],[Hadlat_se1_5mm(1:BHIedy_loc); flipud(Hadlat_se_5mm(1:BHIedy_loc))],deeplavender,'FaceAlpha',0.3,'LineStyle','none')
% patch([Hadyr; flipud(Hadyr)],[hadlat_se1_20ygauss; flipud(hadlat_se_20ygauss)],deeplavender,'FaceAlpha',0.3,'LineStyle','none')
% hold off
% axis tight
% ax65.Box = 'off';
% %ax55.YAxisLocation = 'right';
% ax65.XColor = 'none';
% xlim([1430 2000])
% %ax72.YAxisLocation = 'right';
% ylabel('Hadley Latitude (^oN)')
% %xlabel('Years (CE)')
% ntitle(panel{6},'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% vfill(1797,1820,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% % ax65 = nexttile;
% % dust = rgb('dust');
% % plot(GeoB9501_age,GeoB9501_dustprop,'-','Color',dust,'LineWidth',3)%GeoB9501_dustprop_20ygauss,'-','Color',dust,'LineWidth',3)
% % hold on
% % patch([GeoB9501_age; flipud(GeoB9501_age)],[GeoB9501_dustmax; flipud(GeoB9501_dustmin)],dust,'FaceAlpha',0.3,'LineStyle','none')
% % hold off
% % axis tight
% % ax65.Box = 'off';
% % ax65.XColor = 'none';
% % xlim([1430 2000])
% % ylabel('Aeolian Dust Proportion (%)')
% % ntitle(panel{5},'location','nw')
% % vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% % vfill(1797,1820,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ax75 = nexttile;
% anomaly(YEAR,bhirecon_20ygauss,'thresh',mean(BHIz_ins,'omitnan'))
% axis tight
% ax75.Box = 'off';
% ax75.YAxisLocation = 'right';
% xlim([1430 2000])
% ntitle(panel{7},'location','nw')
% vfill(1645,1715,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% vfill(1797,1820,rgb('light grey blue'),'edgecolor','none','FaceAlpha',0.3,'bottom');
% ylabel('BHI Anomaly')
% xlabel('Year (CE)')
% 
% % Can't wait to probably make more plots...

%
%% Site Coordinates

lat = [36.62, 35.36, 35.68, 35.66, 35.59, 34.34, 34.40, 34.46, 35.52];
lon = [-83.08, -83.93, -83.40, -83.46, -83.92, -84.14, -84.16, -84.18, -83.50];
LITU_coord = [lat; lon]';
QUMO_coord = [35.71; -83.34]';
lat = [35.60, 35.21];
lon = [-83.42, -83.55];
QURU_coord = [lat; lon]';
lat = [30.20, 30.20, 32.39, 29.10, 29.55, 29.44];
lon = [-82.56, -82.56, -84.20, -81.54, -82.37, -82.56];
QUST_coord = [lat; lon]';
lat = [35.30, 37.16, 37.21, 37.17, 35.21, 35.88, 36.22, 35.70, 36.23, 37.17, 34.33, 37.13, 37.17, 37.08, 37.18, 35.49, 35.02];
lon = [-82.38, -80.26, -80.22, -80.28, -83.55, -81.93, -84.08, -84.88, -80.37, -80.17, -87.45, -80.13, -80.48, -83.00, -86.10, -78.41, -83.16];
QUAL_coord = [lat; lon]';
lat = [35.30, 35.37, 32.50, 37.11, 32.50];
lon = [-82.38, -85.26, -84.50, -80.29, -84.50];
QUPR_coord = [lat; lon]';
lat = [35.30, 35.30, 37.11]; lon = [-82.38, -84.19, -80.23];
QSSP_coord = [lat; lon]';
lat = [33.48, 30.28, 32.35, 33.11, 32.05, 35.50, 34.29, 34.95, 33.29, 33.16, 28.18, 30.35, 29.01, 28.55, 32.87, 31.62, 36.45];
lon = [-79.52, -85.53, -81.22, -80.25, -83.30, -76.56, -78.11, -91.22, -79.31, -90.30, -82.05, -88.35, -82.43, -82.16, -89.13, -81.80, -76.62];
TADI_coord = [lat; lon]';
PIPU_coord = [47.15; -80.23]';
lat = [33.49, 30.51, 31.13, 30.20, 32.51, 35.21, 31.01, 35.09, 35.09, 35.40, 35.58, 29.19, 33.55, 27.69, 30.45, 29.84, 29.30, 28.10, 30.72, 32.25, 32.85, 31.07, 35.08, 30.44, 28.15, 27.98];
lon = [-85.42, -83.58, -84.29, -85.15, -84.29, -80.03, -87.15, -79.22, -79.22, -78.54, -78.25, -82.59, -79.14, -81.29, -86.55, -81.95, -83.10, -82.41, -83.25, -85.19, -84.48, -88.98, -79.22, -94.35, -82.08, -81.22];
PIPA_coord = [lat; lon]';
PSSP_coord = [37.11; -80.23]';
lat = [33.11, 31.56, 29.44]; lon = [-80.25, -82.19, -82.57];
QULY_coord = [lat; lon]';
lat = [34.46, 35.45, 35.45, 35.21, 35.17, 36.67, 35.40, 35.26, 35.45, 35.45, 35.08, 34.53, 35.17];
lon = [-84.18, -83.14, -83.14, -83.55, -83.56, -81.67, -83.30, -85.34, -85.57, -83.14, -83.18, -84.39, -82.43];
TSCA_coord = [lat; lon]';
QUBN_coord = [34.04; -86.03]';
QUVE_coord = [35.37; -85.26]';
lat = [32.15, 34.42, 35.49, 35.40, 33.95, 32.25, 35.37]; lon = [-85.12, -82.52, -78.41, -78.54, -83.32, -85.20, -83.56];
PIEC_coord = [lat; lon]';
JUVI_coord = [36.12; -84.00]';
lat = [35.08, 36.24]; lon = [-83.18, -81.33]; 
TSCR_coord = [lat; lon]';
PITA_coord = [35.39; -78.53]';
lat = [30.23, 30.29, 30.02, 30.17, 30.04]; lon = [-85.11, -85.07, -85.12, -85.18, -85.16];
NYOG_coord = [lat; lon]';

%
%% Site Map

% latitude = ncread('p1948.nc','lat');
% longitude = ncread('p1948.nc','lon');
% plon = longitude - 360;

% LITU are circles and yellowish
yellowish = rgb('yellowish');
% Pinus are upward triangles (^)
    % PIPA pine green
    pinegreen = rgb('pine green');
    % PIEC moss green
    mossgreen = rgb('moss green');
    % PIPU dark sea green
    darkseagreen = rgb('dark sea green');
    % PITA dull green
    dullgreen = rgb('dull green');
    % PSSP green grey
    greengrey = rgb('green grey');
% Quercus are squares (s)
    % QUMO is brown grey
    browngrey = rgb('brown grey');
    % QURU is crimson
    crimson = rgb('crimson');
    % QUST is rusty orange
    rustyorange = rgb('rusty orange');
    % QUAL is light grey (I know, right?)
    lightgrey = rgb('light grey');
    % QUPR is chestnut
    chestnut = rgb('chestnut');
    % QSSP is green grey
    % QULY is light beige
    lightbeige = rgb('light beige');
    % QUBN is sandstone
    sandstone = rgb('sandstone');
    % QUVE is black
 % TADI are downward triangles (v) and swamp
 swamp = rgb('swamp');
 % Tsuga are diamonds (d)
    % TSCA is wintergreen 
    wintergreen = rgb('wintergreen');
    % TSCR is cloudy blue
    cloudyblue = rgb('cloudy blue');
 % JUVI are right pointing triangles (>) and blood red
 bloodred = rgb('blood red');
 % NYOG are left pointing triangles (<) and sage green
 sagegreen = rgb('sage green');
 

%
%% Coordinates for used chronologies

QUAL_coord = [35.30 82.38; 37.17 -80.28; 34.33 -87.45; 37.17 -80.48; 37.18 -86.10; 35.49 -78.41; 35.02 -83.16];%[37.17 -80.28; 35.88 -81.93; 37.17 -80.17; 37.13 -80.13; 35.49 -78.41; 35.30 -82.38; 37.17 -80.48];
QULY_coord = [31.56 -82.19; 29.44 -82.57];
TSCA_coord = [34.46 -84.18; 35.26 -85.34; 35.08 -83.18; 34.53 -84.39; 35.17 -82.43];%[35.45 -83.14; 35.08 -83.18; 34.53 -84.39; 35.17 -82.43; 34.46 -84.18; 34.53 -84.39];
TADI_coord = [30.28 -85.53; 33.11 -80.25; 32.05 -83.30; 34.29 -78.11; 34.95 -91.22; 30.35 -88.35; 28.55 -82.16; 32.87 -89.13; 31.62 -81.80];%[34.29 -78.11; 33.16 -90.30; 33.11 -80.25; 28.55 -82.16; 32.87 -89.13; 31.62 -81.80; 33.29 -79.31; 30.35 -88.35; 30.28 -85.53];
PIEC_coord = [35.37 -83.56];%[34.42 -82.52; 33.95 -83.32];
QUMO_coord = [32.25 -85.19];%[32.25 -85.19; 35.71 -83.34];
LITU_coord = [35.59 -83.92; 34.34 -84.14; 34.40 -84.16; 34.46 -84.18];%[35.66 -83.46; 35.59 -83.92; 36.62 -83.08; 35.32 -83.50];
%PIPU_coord = [37.15 -80.23];
PIPA_coord = [27.69 -81.29; 27.98 -81.22];%[33.49 -85.42; 30.51 -83.58; 35.21 -80.03; 35.09 -79.22; 35.09 -79.22; 33.55 -79.14; 30.72 -83.25; 32.51 -84.29; 31.07 -88.98; 29.84 -81.95];
QUPR_coord = [35.37 -85.26; 32.50 -84.50; 37.11 -80.29]%[32.50 -84.50; 37.11 -80.29];
%PITA_coord = [35.39 -78.53];
QUSP_coord = [35.30 -82.38];%[35.30 -82.38; 34.28 -84.19];
QUST_coord = [32.39 -84.20; 29.55 -82.37];%[29.55 -82.37];
%NYOG_coord = [30.17 -85.18];
%QUVE_coord = [35.37 -85.26]; 
QUBN_coord = [34.04 -86.03];
JUVI_coord = [36.12 -84.00];

%
%% Just shortening the run time


figure
%usbasemap
geobasemap darkwater
ax = gca;
%gx = geoaxes;
%ax.Color = rgb('pale aqua');
%xlim([-97 -70])
%ylim([24 40])
geolimits([24 38],[-92 -75]);
%yellowish = rgb('yellowish');

hold on
geoscatter(JUVI_coord(:,1),JUVI_coord(:,2),300,'p','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('blood red'));
hold off
%yellowish = rgb('yellowish');
hold on
geoscatter(LITU_coord(:,1),LITU_coord(:,2),300,'o','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('yellowish'));
hold off

% hold on
% geoscatter(NYOG_coord(:,1),NYOG_coord(:,2),100,'h','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('clay')); %rgb('sage green'));
% hold off

hold on
geoscatter(PIEC_coord(:,1),PIEC_coord(:,2),300,'^','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('mossy green'));
hold off

hold on
geoscatter(PIPA_coord(:,1),PIPA_coord(:,2),300,'^','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('pine'));
hold off

% hold on
% geoscatter(PIPU_coord(:,1),PIPU_coord(:,2),100,'^','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('jade')); %rgb('dark sea green'));
% hold off

% hold on
% geoscatter(PITA_coord(:,1),PITA_coord(:,2),100,'^','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('emerald'));%rgb('dull green'));
% hold off

% hold on
% geoscatter(PSSP_coord(:,1),PSSP_coord(:,2),100,'^','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('green grey'));
% hold off

hold on
geoscatter(QUAL_coord(:,1),QUAL_coord(:,2),300,'s','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('carolina blue'));%rgb('off white'));%light grey'));
hold off

hold on
geoscatter(QUBN_coord(:,1),QUBN_coord(:,2),300,'s','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('cerulean'));%rgb('sandstone'));
hold off

hold on
geoscatter(QULY_coord(:,1),QULY_coord(:,2),300,'s','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('blue')); %rgb('light beige'));
hold off

hold on
geoscatter(QUMO_coord(:,1),QUMO_coord(:,2),300,'s','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('water blue'));%rgb('dark aquamarine'));%brown grey'));
hold off

hold on
geoscatter(QUPR_coord(:,1),QUPR_coord(:,2),300,'s','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('electric blue'));%rgb('chestnut'));
hold off

% hold on
% geoscatter(QURU_coord(:,1),QURU_coord(:,2),100,'s','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('crimson'));
% hold off

hold on
geoscatter(QUST_coord(:,1),QUST_coord(:,2),300,'s','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('medium blue'));%rgb('rusty orange'));
hold off

% hold on
% geoscatter(QUVE_coord(:,1),QUVE_coord(:,2),100,'s','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('deep sky blue'));%'k');
% hold off

hold on
geoscatter(QUSP_coord(:,1),QUSP_coord(:,2),300,'s','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('neon blue'));%rgb('blurple'));%green grey'));
hold off

hold on
geoscatter(TADI_coord(:,1),TADI_coord(:,2),300,'v','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('swamp'));
hold off

hold on
geoscatter(TSCA_coord(:,1),TSCA_coord(:,2),300,'d','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('orchid'));
hold off

% hold on
% geoscatter(PIPU_coord(:,1),PIPU_coord(:,2),100,'^','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('jade')); %rgb('dark sea green'));
% hold off
% hold on
% geoscatter(TSCR_coord(:,1),TSCR_coord(:,2),100,'d','filled','MarkerEdgeColor','k','MarkerFaceColor',rgb('cloudy blue'));
% hold off

legend('{\it J. virginiana}','{\it L. tulipifera}','{\it P. echinata}','{\it P. palustris}','{\it Q. alba}','{\it Q. boyntonii}','{\it Q. lyrata}','{\it Q. montana}','{\it Q. prinus}','{\it Q. stellata}','{\it Q. ssp}','{\it T. distichum}','{\it T. canadensis}','Location','southeast')

%legend('{\it L. tulipifera}','{\it P. palustris}','{\it P. echinata}','{\it P. pungens}','{\it P. taeda}','{\it P. ssp}','{\it Q. montana}','{\it Q. rubra}','{\it Q. stellata}','{\it Q. alba}','{\it Q. prinus}','{\it Q. ssp}','{\it Q. lyrata}','{\it Q. boyntonii}','{\it Q. velutina}','{\it T. distichum}','{\it T. canadensis}','{\it T. caroliniana}','{\it J. virginiana}','{\it N. ogeche}','Location','southeast')
%ax.Color = rgb('very light blue');%'eggshell');
%ax.AmbientLightColor = rgb('very light blue');
%usbasemap
%ax.Colormap = rgb('eggshell');
%colormap(rgb('eggshell'))
%ay = gca;
%ay.Color = rgb('eggshell');
%xlim([-95 -70])
%ylim([24 38])
%geolimits([24 38],[-92 -75]);

ax.Box = 'on';
%ax.Color = rgb('pale aqua');
%axis tight

% I don't really need anything beyond this point.

%
%% Wavelet Analysis

[twave,tper,tscale,tcoi,tsig95] = wt([RECON],'Mother','morlet','MakeFigure',1);


%
%% Autocorrelation

BHI_acf = acf(RECON,25);


%
%% PDSI Connections

NADAtime = ncread('NADAv2-2008.nc','time'); NADAlat = ncread('NADAv2-2008.nc','lat'); NADAlon = ncread('NADAv2-2008.nc','lon'); NADAPDSI = ncread('NADAv2-2008.nc','PDSI');
SE_NADA_PDSI = NADAPDSI(18:26,11:16,1:577);
SE_NADA_PDSI_r = nan(9,6);
SE_NADA_PDSI_p = nan(9,6);
for i = 1:9
for j = 1:6
[rnada,pnada] = corr(flipud(RECON(1:577)),squeeze(SE_NADA_PDSI(i,j,:)));
SE_NADA_PDSI_r(i,j) = rnada;
SE_NADA_PDSI_p(i,j) = pnada;
end
end

SE_NADA_lat_rep = repmat(NADAlat(11:16),[1 9])';
 SE_NADA_PDSI_mean = nan(1,577)';
for i = 1:577
SE_NADA_PDSI_mean(i) = wmean(squeeze(SE_NADA_PDSI(:,:,i)),cosd(SE_NADA_lat_rep),'all','omitnan');
end
[rnadamean,pnadamean] = corr(flipud(RECON(1:577)),SE_NADA_PDSI_mean,'tail','right');

figure
yyaxis left
plot(YEAR,RECON,'k-')
yyaxis right
plot(YEAR(1:577),flipud(SE_NADA_PDSI_mean),'r-')
xlim([1430 2000])

%
%% Standardized Precipitation Index Connections


SPI_lat = ncread('NASPA_WARM_SPI.nc','lat'); 
SPI_time = ncread('NASPA_WARM_SPI.nc','time'); 
SPI_lon = ncread('NASPA_WARM_SPI.nc','lon'); 
SPI_SPI = ncread('NASPA_WARM_SPI.nc','SPI');
SPI_SE_lat = SPI_lat(26:49); SPI_SE_lon = SPI_lon(83:115); SPI_SE_time = SPI_time(1431:2012); 
SPI_SE_SPI = SPI_SPI(1431:2012,26:49,83:115);

SPI_SE_lat_rep = repmat(SPI_SE_lat,[1 33]);
SPI_SE_SPI_mean = nan(1,582);
for i = 1:582
    SPI_SE_SPI_mean(i) = wmean(squeeze(SPI_SE_SPI(i,:,:)),cosd(SPI_SE_lat_rep),'all','omitnan');
end
[rspimean,pspimean] = corr(RECON,SPI_SE_SPI_mean','tail','right');

% Total precipitation (from SPI?)
SPI_total_precip = ncread('NASPA_WARM_TOTALPRECIP.nc','precip');

SPI_SE_total_precip = SPI_total_precip(1431:2012,26:49,83:115);
SPI_SE_total_precip_mean = nan(1,582)';
for i = 1:582
    SPI_SE_total_precip_mean(i) = wmean(squeeze(SPI_SE_total_precip(i,:,:)),cosd(SPI_SE_lat_rep),'all','omitnan');
end
[rspitpmean,pspitpmean] = corr(RECON,SPI_SE_total_precip_mean,'tail','right');

figure
subplot(3,2,[1,2])
% panel = {'A','B','C','D','E','F'};
% %midnight = rgb('midnight');
% %ultramar = rgb('ultramarine');
% BHIz_ins_10mm = movmean(BHIz_ins,10,'EndPoints','fill');
% BHIz_ins_5mm = movmean(BHIz_ins,5,'EndPoints','fill');
% figure
% %subplot(6,2,[1,2])
% ax(1) = subsubplot(6,1,1,'vpad',0.01);
% %ax = gca;
waterblue = rgb('water blue');
plot(YEAR,RECON,'-','Color',waterblue,'LineWidth',2)
% hold on
% plot(YEAR,BHIz_recon_10movmean,'-','Color',midnight,'LineWidth',3)
% hold off
hold on
patch([YEAR; flipud(YEAR)],[PREUP_95; flipud(PRELO_05)],waterblue,'FaceAlpha',0.3,'LineStyle','none')
hold off
axis tight
ylabel('BHI')
% ax(1).XAxisLocation = 'top';
% ax(1).XColor = 'k';
% ntitle(panel{1},'location','nw')
subplot(3,2,[3,4])
plot(YEAR(1:577),flipud(SE_NADA_PDSI_mean),'r-','LineWidth',2)
axis tight
ylabel('PDSI')
subplot(3,2,[5,6])
ax = gca;
yyaxis left
plot(YEAR,SPI_SE_SPI_mean,'k-','LineWidth',2)
ylabel('SPI')
ax.YColor = 'k';
yyaxis right
plot(YEAR,SPI_SE_total_precip_mean,'g-','LineWidth',2)
ylabel('Total Precipitation (mm)')
ax.YColor = 'k';
xlabel('Year (CE)')
axis tight

%
%% Stats for pre-industrial vs industrial periods

% Added 7/26/2022
preindustrial_recon = nan(length(RECON),1);
preindustrial_recon(1:find(YEAR == 1849,1,'first'),1) = RECON(1:find(YEAR == 1849,1,'first'),1);
industrial_recon = nan(length(RECON),1);
industrial_recon(find(YEAR == 1850,1,'first'):end,1) = RECON(find(YEAR == 1850,1,'first'):end,1);
recon_comp = [preindustrial_recon,industrial_recon];

[comp_p,comp_stats] = vartestn(recon_comp,'TestType','Bartlett')

[panova1_comp,tblanova1_comp,statsanova1_comp] = anova1(recon_comp)

[p_med_comp,h_med_comp,stats_med_comp] = ranksum(preindustrial_recon,industrial_recon,'alpha',0.05,'method','approximate','tail','both')

p90 = prctile(RECON,90);
p95 = prctile(RECON,95);
p99 = prctile(RECON,99);
p10 = prctile(RECON,10);
p05 = prctile(RECON,5);
p01 = prctile(RECON,1);

BHI_90 = RECON;
BHI_90(BHI_90 < p90)=NaN;
BHI_95 = RECON;
BHI_95(BHI_95 < p95) = NaN;
BHI_99 = RECON;
BHI_99(BHI_99 < p99) = NaN;
BHI_10 = RECON;
BHI_10(BHI_10 > p10) = NaN;
BHI_05 = RECON;
BHI_05(BHI_05 > p05) = NaN;
BHI_01 = RECON;
BHI_01(BHI_01 > p01) = NaN;