% Figures for MEPS paper October 2023
% Map of West Antarctic Peninsula
% Map info thanks to Schuyler Nardelli & Jack Conroy, March 2020
% Updated by Jessie Turner June 2021 through October 2023
% Figure 5 in new submission to MEPS

clear; close all hidden;
% Load data, saved from script "calculate_bloomIndices_WAP.m"
datafolder = '/Users/jsturner/satellite_data/matlab_WAP/daily_CMEMS/';
load([datafolder 'StartDates_CMEMS_WAP_v4.mat']);
% Load map information
mapfolder = '/Users/jsturner/satellite_data/matlab_WAP/WAP_map_SchuylerJack/';
load([mapfolder 'bathy_and_coast.mat'],'coastlon','coastlat','shlflon','shlflat','lon','lat','depth'); % load Jack+Schuyler's coastline:
load([mapfolder 'WAP_regions_2022-07-15.mat']);
load([mapfolder 'continental_shelf_WAP.mat']);
% Fix peak date to be 365 days and more rather than 1-180.
tmax(tmax<180) = tmax(tmax<180)+365;
reg_str_list = {'spf','miz','nos','mis','sos'};

%%  Average over decades for bloom start date and peak date
% ind_0010 = find(yrs >= 2001 & yrs < 2011); % Jan 2001 to Dec 2010
% ind_1020 = find(yrs >= 2011 & yrs < 2021); % Jan 2011 to Dec 2020
ind_0010 = find(yrs >= 2001 & yrs < 2012); % Jan 2001 to Jan 2011
ind_1020 = find(yrs >= 2012 & yrs < 2023); % Jan 2012 to Jan 2022

bi_1 = squeeze(bi_thres(:,:,ind_0010));
bi_avg_thres_1 = nanmean(bi_1,3);
bi_2 = squeeze(bi_thres(:,:,ind_1020));
bi_avg_thres_2 = nanmean(bi_2,3);

tmax_1 = squeeze(tmax(:,:,ind_0010));
tmax_avg_1 = nanmean(tmax_1,3);
tmax_2 = squeeze(tmax(:,:,ind_1020));
tmax_avg_2 = nanmean(tmax_2,3);

% Difference between decades
bi_decade_diff = bi_avg_thres_2 - bi_avg_thres_1;
tmax_decade_diff = tmax_avg_2 - tmax_avg_1;

% Set up figure
figure(1)
set(gcf,'position',[8         200        1000         600]);
clf

%% Start date
% Set up ticks for day of year
dates_ticks = [244 274 305 335 366 397];
datestr_bi = list_doy_strings(dates_ticks);

ax(1) = subplot(2,3,1);
m_proj('lambert','long',[-80 -55],'lat',[-70.05 -59.95]); hold on
m_pcolor(lon_chl,lat_chl,bi_avg_thres_1); shading flat;
colormap(ax(1),turbo(30));
hold on
m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,...
    'xticklabels',[],'yticklabels',[],'fontsize',16,'linewidth',1,...
        'fontname','Arial','backcolor','none');caxis([dates_ticks(1) dates_ticks(end)]);
caxis([dates_ticks(1) dates_ticks(end)]);
% Add coast
coastcolor = [.8 .8 .8];
for i=1:length(coastlon)
    p_coast = m_patch(coastlon{i},coastlat{i},coastcolor,'LineWidth',1);
end; hold on;
% Add ice shelf
shelfcolor = [1 1 1];
for i=1:length(shlflon)
    p_shelf = m_patch(shlflon{i},shlflat{i},shelfcolor,'LineWidth',1);
end
m_plot(contshlf_lon,contshlf_lat,'linestyle',':','color','k','linewidth',2);
title('2001-2011','fontsize',16,'fontweight','normal');
pos = get(gca,'position');
set(gca,'position',[pos(1)-0.05,pos(2)+0.03,pos(3)+0.055,pos(4)]);


ax(2) = subplot(2,3,2);
m_proj('lambert','long',[-80 -55],'lat',[-70.05 -59.95]); hold on
m_pcolor(lon_chl,lat_chl,bi_avg_thres_2); shading flat;
colormap(ax(2),turbo(30));
hold on
m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,...
    'xticklabels',[],'yticklabels',[],'fontsize',16,'linewidth',1,...
        'fontname','Arial','backcolor','none');caxis([dates_ticks(1) dates_ticks(end)]);
caxis([dates_ticks(1) dates_ticks(end)]);
cb2 = colorbar;
cb2.Location = 'southoutside';
cb2.Label.String = 'Start Date';
cb2.LineWidth = 1.5;
cb2.Ticks = dates_ticks;
cb2.TickLabels = {'S','O','N','D','J','F'};
cb2.FontSize = 16;
cb2.Position = [0.25    0.57    0.2    0.03];
% Add coast
coastcolor = [.8 .8 .8];
for i=1:length(coastlon)
    p_coast = m_patch(coastlon{i},coastlat{i},coastcolor,'LineWidth',1);
end; hold on;
% Add ice shelf
shelfcolor = [1 1 1];
for i=1:length(shlflon)
    p_shelf = m_patch(shlflon{i},shlflat{i},shelfcolor,'LineWidth',1);
end
m_plot(contshlf_lon,contshlf_lat,'linestyle',':','color','k','linewidth',2);
title('2012-2022','fontsize',16,'fontweight','normal');
pos = get(gca,'position');
set(gca,'position',[pos(1)-0.05,pos(2)+0.03,pos(3)+0.055,pos(4)]);


ax(3) = subplot(2,3,3); cla;
m_proj('lambert','long',[-80 -55],'lat',[-70.05 -59.95]); hold on
m_pcolor(lon_chl,lat_chl,bi_decade_diff); shading flat;
colormap(ax(3),flipud(cbrewer('div','RdBu',30)));
hold on
m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,...
    'xticklabels',[],'yticklabels',[],'fontsize',16,'linewidth',1,...
        'fontname','Arial','backcolor','none');caxis([dates_ticks(1) dates_ticks(end)]);
caxis([-40 40]);
cb3 = colorbar;
cb3.Location = 'southoutside';
cb3.Label.String = 'Difference 2012-2022 minus 2001-2011 (days)';
cb3.LineWidth = 1.5;
cb3.FontSize = 16;
cb3.Position = [0.675    0.57    0.2    0.03];
% Add coast
coastcolor = [.8 .8 .8];
for i=1:length(coastlon)
    p_coast = m_patch(coastlon{i},coastlat{i},coastcolor,'LineWidth',1);
end; hold on;
% Add ice shelf
shelfcolor = [1 1 1];
for i=1:length(shlflon)
    p_shelf = m_patch(shlflon{i},shlflat{i},shelfcolor,'LineWidth',1);
end
m_plot(contshlf_lon,contshlf_lat,'linestyle',':','color','k','linewidth',2);
% title({'Difference 2010-2020 minus 2000-2010';''},'fontsize',16,'fontweight','normal');
title(' ');
pos = get(gca,'position');
set(gca,'position',[pos(1)-0.05,pos(2)+0.03,pos(3)+0.055,pos(4)]);

%% Peak date
% Set up ticks for day of year
dates_ticks2 = [274 305 335 366 397 425 456];
datestr_tmax = list_doy_strings(dates_ticks2);

ax(4) = subplot(2,3,4);
m_proj('lambert','long',[-80 -55],'lat',[-70.05 -59.95]); hold on
m_pcolor(lon_chl,lat_chl,tmax_avg_1); shading flat;
colormap(ax(4),turbo(30));
hold on
m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,...
    'xticklabels',[],'yticklabels',[],'fontsize',16,'linewidth',1,...
        'fontname','Arial','backcolor','none');caxis([dates_ticks2(1) dates_ticks2(end)]);
% Add coast
coastcolor = [.8 .8 .8];
for i=1:length(coastlon)
    p_coast = m_patch(coastlon{i},coastlat{i},coastcolor,'LineWidth',1);
end; hold on;
% Add ice shelf
shelfcolor = [1 1 1];
for i=1:length(shlflon)
    p_shelf = m_patch(shlflon{i},shlflat{i},shelfcolor,'LineWidth',1);
end
m_plot(contshlf_lon,contshlf_lat,'linestyle',':','color','k','linewidth',2);
title('2001-2011','fontsize',16,'fontweight','normal');
pos = get(gca,'position');
set(gca,'position',[pos(1)-0.05,pos(2)+0.01,pos(3)+0.055,pos(4)]);


ax(5) = subplot(2,3,5);
m_proj('lambert','long',[-80 -55],'lat',[-70.05 -59.95]); hold on
m_pcolor(lon_chl,lat_chl,tmax_avg_2); shading flat;
colormap(ax(5),turbo(30));
hold on
m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,...
    'xticklabels',[],'yticklabels',[],'fontsize',16,'linewidth',1,...
        'fontname','Arial','backcolor','none');caxis([dates_ticks2(1) dates_ticks2(end)]);
cb4 = colorbar;
cb4.Location = 'southoutside';
cb4.Label.String = 'Peak Date';
cb4.LineWidth = 1.5;
cb4.Ticks = dates_ticks2;
cb4.TickLabels = {'O','N','D','J','F','M','A'};
cb4.FontSize = 16;
cb4.Position = [0.25    0.075    0.2    0.03];
% Add coast
coastcolor = [.8 .8 .8];
for i=1:length(coastlon)
    p_coast = m_patch(coastlon{i},coastlat{i},coastcolor,'LineWidth',1);
end; hold on;
% Add ice shelf
shelfcolor = [1 1 1];
for i=1:length(shlflon)
    p_shelf = m_patch(shlflon{i},shlflat{i},shelfcolor,'LineWidth',1);
end
m_plot(contshlf_lon,contshlf_lat,'linestyle',':','color','k','linewidth',2);
title('2012-2022','fontsize',16,'fontweight','normal');
pos = get(gca,'position');
set(gca,'position',[pos(1)-0.05,pos(2)+0.01,pos(3)+0.055,pos(4)]);


ax(6) = subplot(2,3,6); cla;
m_proj('lambert','long',[-80 -55],'lat',[-70.05 -59.95]); hold on
m_pcolor(lon_chl,lat_chl,tmax_decade_diff); shading flat;
colormap(ax(6),flipud(cbrewer('div','RdBu',30)));
hold on
m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,...
    'xticklabels',[],'yticklabels',[],'fontsize',16,'linewidth',1,...
        'fontname','Arial','backcolor','none');
caxis([-50 50]);
cb5 = colorbar;
cb5.Location = 'southoutside';
cb5.Label.String = 'Difference 2012-2022 minus 2001-2011 (days)';
cb5.LineWidth = 1.5;
cb5.FontSize = 16;
cb5.Position = [0.675    0.075    0.2    0.03];
% Add coast
coastcolor = [.8 .8 .8];
for i=1:length(coastlon)
    p_coast = m_patch(coastlon{i},coastlat{i},coastcolor,'LineWidth',1);
end; hold on;
% Add ice shelf
shelfcolor = [1 1 1];
for i=1:length(shlflon)
    p_shelf = m_patch(shlflon{i},shlflat{i},shelfcolor,'LineWidth',1);
end
m_plot(contshlf_lon,contshlf_lat,'linestyle',':','color','k','linewidth',2);
% title({'Difference 2010-2020 minus 2000-2010';''},'fontsize',16,'fontweight','normal');
title(' ');
pos = get(gca,'position');
set(gca,'position',[pos(1)-0.05,pos(2)+0.01,pos(3)+0.055,pos(4)]);


%% Save figure

sublabel(ax,2,14,'fontsize',18);

print(gcf,'Fig_5_Decadal_Maps_StartDate_PeakDate_11yrs.png','-dpng','-r200');
print(gcf,'Fig_5_Decadal_Maps_StartDate_PeakDate_11yrs.eps','-depsc');




