%%%% FIGURE 3 for paper
% Explore CMEMS data - GlobColour
% Data source:
% https://resources.marine.copernicus.eu/product-detail/OCEANCOLOUR_GLO_CHL_L4_REP_OBSERVATIONS_009_082/DATA-ACCESS
% As of October 2023

clear; close all hidden; format shortg;
% 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([datafolder 'Mean_Chl_20231026.mat']);
% Load map information
mapfolder = '/Users/jsturner/satellite_data/matlab_WAP/WAP_map_SchuylerJack/';
load([mapfolder 'bathy_and_coast.mat'],'coastlon','coastlat','shlflon','shlflat'); % load Jack+Schuyler's coastline:
load([mapfolder 'WAP_regions_2022-07-15.mat']);
load([mapfolder 'continental_shelf_WAP.mat']);

%% Average over all years for all bloom indices
% Long-term means of all indices
bi_avg_thres = nanmean(bi_thres,3);
% Fix Tmax to be 365 days and more rather than 1-180.
tmax(tmax<180) = tmax(tmax<180)+365;
tmax_avg = nanmean(tmax,3);

%% Plot each bloom index long-term average start date (map)
close all hidden;
% Set up ticks for day of year
dates_ticks = [244 274 305 335 366 397 425 456];
figure('position',[291   229   692   334]);
clf

%%%%%% a. Start Date (Threshold method)
haxes(1) = subplot(1,2,1);
m_proj('lambert','long',[-80 -55],'lat',[-70 -60]); hold on
m_pcolor(lon_chl,lat_chl,bi_avg_thres); shading flat;
m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,'fontsize',14,'linewidth',1,'fontname','Arial','backcolor','none','xticklabels',[],'yticklabels',[]);
colormap(turbo(30));
caxis([dates_ticks(1) dates_ticks(end-1)]);
cb = colorbar;
cb.Label.String = 'Start Date';
cb.LineWidth = 1.5;
cb.Ticks = dates_ticks;
cb.TickLabels = {'S','O','N','D','J','F','M'};
cb.Location = 'southoutside';
cb.FontSize = 14;
cb.Label.FontSize = 16;
cb.Position = [0.15 0.15 0.295 0.047904];
hold on
% 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
% Add continental shelf break
m_plot(contshlf_lon,contshlf_lat,'linestyle',':','color','k','linewidth',1.5);
pos = get(gca,'position');
set(gca,'position',[pos(1)-0.04,pos(2)-0.04,pos(3)+0.08,pos(4)+0.08]);

%%%%%% b. Timing of maximum Chlorophyll "Peak date"
haxes(2) = subplot(1,2,2);
m_proj('lambert','long',[-80 -55],'lat',[-70 -60]); hold on
m_pcolor(lon_chl,lat_chl,tmax_avg); shading flat;
m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,'fontsize',14,'linewidth',1,'fontname','Arial','backcolor','none','xticklabels',[],'yticklabels',[]);
colormap(turbo(30));
caxis([dates_ticks(3) dates_ticks(end)]);
cb2 = colorbar;
cb2.Label.String = 'Peak Date';
cb2.LineWidth = 1.5;
cb2.Ticks = dates_ticks(3:end);
cb2.TickLabels = {'N','D','J','F','M','A'};
cb2.Location = 'southoutside';
cb2.FontSize = 14;
cb2.Label.FontSize = 16;
cb2.Position = [0.59 0.15 0.295 0.047904];
hold on
% 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
% Add continental shelf break
m_plot(contshlf_lon,contshlf_lat,'linestyle',':','color','k','linewidth',1.5);
pos = get(gca,'position');
set(gca,'position',[pos(1)-0.04,pos(2)-0.04,pos(3)+0.08,pos(4)+0.08]);

sublabel(haxes,1,14);

% SAVE FIGURE
print(gcf,'Fig_3_Avg_Bloom_Indices_twoSubplots.png','-dpng','-r300');
print(gcf,'Fig_3_Avg_Bloom_Indices_twoSubplots.eps','-depsc');
