% Explore CMEMS data - GlobColour
% Sources
% https://resources.marine.copernicus.eu/product-detail/OCEANCOLOUR_GLO_CHL_L4_REP_OBSERVATIONS_009_082/DATA-ACCESS
% Jessie Turner January 2022
% Updated to plot a time series of L4 chl from daily interpolated product
% Updated May 11, 2022
% To compare 5 bloom start date indices, Threshold is only one.
% Compare trends in the different bloom start dates among indices.
% Updated January 2024 to use Theil-Sen Slope

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 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']);

%% Calculate trend over time for all bloom indices
years = [1998:2022]';
% Fix Tmax to be 365 days and more rather than 1-180.
tmax(tmax<180) = tmax(tmax<180)+365;

% Start Date
% Old method: Calculate trend ignoring NaNs, calculate significance of trend from Mann-Kendall test.
% [trendmap_bi,tau_bi] = trend(bi_thres,years,'omitnan','Type','Kendall');
% [~,pval_bi] = mann_kendall(bi_thres);
% compute trend over years at each x.y point Theil-Sen
disp('Calculating start date trend at each x,y point (may be slow)...');
bi_data = reshape(bi_thres,[size(bi_thres,1).*size(bi_thres,2),length(years)]);
trendmap_bi_raw = nan(length(bi_data),1);
pval_bi_raw = nan(length(bi_data),1);
for ii = 1:size(bi_data,1)
    data = [years,bi_data(ii,:)'];
    [m,b] = TheilSen(data);
    line0 = years.* m + b;
    % [~,pval0] = corr(data); % test whether the time axis and the data are correlated
    [~,pval0] = corr([line0 bi_data(ii,:)']); % test whether the line and the data are correlated
    trendmap_bi_raw(ii) = m;
    pval_bi_raw(ii) = pval0(2);
end
trendmap_bi = reshape(trendmap_bi_raw,[size(bi_thres,1),size(bi_thres,2)]);
pval_bi = reshape(pval_bi_raw,[size(bi_thres,1),size(bi_thres,2)]);
sigpts_bi = pval_bi < 0.05;

% Peak Date
% Old method: Calculate trend ignoring NaNs, calculate significance of trend from Mann-Kendall test.
% [trendmap_tmax,tau_tmax] = trend(tmax,years,'omitnan','Type','Kendall');
% [~,pval_tmax] = mann_kendall(tmax);
% compute trend over years at each x.y point Theil-Sen
disp('Calculating peak date trend at each x,y point (may be slow)...');
tmax_data = reshape(tmax,[size(tmax,1).*size(tmax,2),length(years)]);
trendmap_tmax_raw = nan(length(tmax_data),1);
pval_tmax_raw = nan(length(tmax_data),1);
for ii = 1:size(tmax_data,1)
    data = [years,tmax_data(ii,:)'];
    [m,b] = TheilSen(data);
    line0 = years.* m + b;
    % [~,pval0] = corr(data); % test whether the time axis and the data are correlated
    [~,pval0] = corr([line0 tmax_data(ii,:)']); % test whether the line and the data are correlated
    trendmap_tmax_raw(ii) = m;
    pval_tmax_raw(ii) = pval0(2);
end
trendmap_tmax = reshape(trendmap_tmax_raw,[size(tmax,1),size(tmax,2)]);
pval_tmax = reshape(pval_tmax_raw,[size(tmax,1),size(tmax,2)]);
sigpts_tmax = pval_tmax < 0.05;


%% Plot each bloom index trend (map)
close all hidden;
figure('position',[291   229   692   334]);
clf

%%%%%% a. Start Date (Threshold)
haxes(1) = subplot(1,2,1);
m_proj('lambert','long',[-80 -55],'lat',[-70 -60]); hold on
m_pcolor(lon_chl,lat_chl,trendmap_bi); 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(flipud(cbrewer('div','RdBu',31)));
caxis([-5 5]);
hold on
[xs,ys] = m_ll2xy(lon_chl,lat_chl);
stipple(xs,ys,sigpts_bi,'markersize',8);
cb = colorbar;
cb.Label.String = {'Trend in Start Date (days yr^{-1})';'::: p < 0.05'};
cb.LineWidth = 1;
cb.Location = 'southoutside';
cb.FontSize = 14;
cb.Label.FontSize = 16;
cb.Position = [0.15 0.21 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),pos(3)+0.08,pos(4)+0.04]);


% Calculate percent of the map experiencing later timing
later_bi = trendmap_bi(trendmap_bi>0);
perc_later = numel(later_bi) ./ numel(trendmap_bi) .* 100;

%%%%% b. Timing of Maximum
haxes(2) = subplot(1,2,2);
m_proj('lambert','long',[-80 -55],'lat',[-70 -60]); hold on
m_pcolor(lon_chl,lat_chl,trendmap_tmax); 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(flipud(cbrewer('div','RdBu',31)));
caxis([-10 10]);
hold on
[xs,ys] = m_ll2xy(lon_chl,lat_chl);
stipple(xs,ys,sigpts_tmax,'markersize',8);
cb2 = colorbar;
cb2.Label.String = {'Trend in Peak Date (days yr^{-1})';'::: p < 0.05'};
cb2.LineWidth = 1;
cb2.Location = 'southoutside';
cb2.FontSize = 14;
cb2.Label.FontSize = 16;
cb2.Position = [0.59 0.21 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),pos(3)+0.08,pos(4)+0.04]);


sublabel(haxes,1,14);
% 
%% SAVE FIGURE
% 
print(gcf,'Fig_8_Trends_Bloom_Indices_twoSubplots.png','-dpng','-r200');
print(gcf,'Fig_8_Trends_Bloom_Indices_twoSubplots.eps','-depsc');
