% Look at satellite ocean color data for West Antarctic Peninsula
% For NASA IDS WAP Project
% Coordinates of ocean color data: W: -80, N: -60, E:-55, S:-70
% Original code from June-July 2021 for MODIS-Aqua only. 
% Updated by Jessie Turner October 2023
clear; close all hidden;

%% Load data
% Load monthly chl file 
disp(' ');
load('/Users/jsturner/satellite_data/matlab_WAP/means_CMEMS/Month_Means_v3.mat');
disp('File loaded.');
% Load info for map
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 'continental_shelf_WAP.mat']);
%% Clean up data
 % Remove pixels that have < 30% of monthly averages
 % Calculate how many points are present for each pixel
time_len = size(md_chl,3); % number of timesteps 
perc_cov = sum(isfinite(md_chl),3)./time_len;
% Trim to only analyze pixels that have >60% of the record.
chl_data = reshape(md_chl,[(size(perc_cov,1).*size(perc_cov,2)),time_len]);
chl_data((perc_cov <= 0.3),:) = nan;
perc_cov_data = reshape(perc_cov,(size(perc_cov,1).*size(perc_cov,2)),1);
perc_cov_data(perc_cov <= 0.3) = nan;
md_chl = reshape(chl_data,[size(perc_cov,1),size(perc_cov,2),time_len]);
 
%% Plot all month trends
close all hidden;
mostr = {'September','October','November','December',...
    'January','February','March','April',};
date_columns = datevec(monthly_date);
yy = date_columns(:,1);
mm = date_columns(:,2);
month_ind = [9 10 11 12 1 2 3 4];

% Calculate long-term mean of each month for relative trend
monthclim_chl = [];
for k = 1:12
    ind_thismo = find(mm==k);
    chl_thismo = squeeze(md_chl(:,:,ind_thismo));
    % how many points are included (data availability?)
    % Calculate how many points are present for each pixel
    time_len = size(chl_thismo,3); % number of timesteps
    perc_cov = sum(isfinite(chl_thismo),3)./time_len;
    % Trim to only analyze pixels that have >60% of the years for that month.
    chl_data = reshape(chl_thismo,[(size(perc_cov,1).*size(perc_cov,2)),time_len]);
    chl_data((perc_cov <= 0.75),:) = nan;
    perc_cov_data = reshape(perc_cov,(size(perc_cov,1).*size(perc_cov,2)),1);
    perc_cov_data(perc_cov <= 0.75) = nan;
    chl_trimmed = reshape(chl_data,[size(perc_cov,1),size(perc_cov,2),time_len]);
    mean_thismo = nanmean(chl_trimmed,3);
    monthclim_chl = cat(3,monthclim_chl,mean_thismo);
end

%% Calculate Trends
trendmap_rel_stack = nan(size(monthclim_chl,1),size(monthclim_chl,2),length(month_ind));
pval_stack = nan(size(monthclim_chl,1),size(monthclim_chl,2),length(month_ind));
for i = 1:length(month_ind)
    disp(char(mostr(i)));
    m1 = month_ind(i);
    ind_thismonth = find(mm == m1);
    chl_thismonth = md_chl(:,:,ind_thismonth);
    chlmean_thismonth = monthclim_chl(:,:,m1);
    
    % Calculate how many points are present for each pixel
    time_len = size(chl_thismonth,3); % number of timesteps
    perc_cov = sum(isfinite(chl_thismonth),3)./time_len;
    % Trim to only analyze pixels that have >75% of the years for that month.
    chl_data = reshape(chl_thismonth,[(size(perc_cov,1).*size(perc_cov,2)),time_len]);
    chl_data((perc_cov <= 0.75),:) = nan;
    perc_cov_data = reshape(perc_cov,(size(perc_cov,1).*size(perc_cov,2)),1);
    perc_cov_data(perc_cov <= 0.75) = nan;
    chl_trimmed = reshape(chl_data,[size(perc_cov,1),size(perc_cov,2),time_len]);
    
    
    t2 = (1:1:size(chl_thismonth,3))';
%     [trendmap_thismonth,~] = trend(chl_trimmed,t2,'omitnan','Type','Kendall');
%     [~,pval] = mann_kendall(chl_trimmed);
    % compute trend over years at each x.y point
    disp('Calculating trend at each x,y point (may be slow)...');
    trendmap_thismonth_raw = nan(length(chl_data),1);
    pval_thismonth_raw = nan(length(chl_data),1);
    for ii = 1:size(chl_data,1)
        data = [t2,chl_data(ii,:)'];
        [m,b] = TheilSen(data);
        line0 = t2.* m + b;
        % [~,pval0] = corr(data); % test whether the time axis and the data are correlated
        [~,pval0] = corr([line0 chl_data(ii,:)']); % test whether the line and the data are correlated
        trendmap_thismonth_raw(ii) = m;
        pval_thismonth_raw(ii) = pval0(2);
    end
    trendmap_thismonth = reshape(trendmap_thismonth_raw,[size(perc_cov,1),size(perc_cov,2)]);
    pval = reshape(pval_thismonth_raw,[size(perc_cov,1),size(perc_cov,2)]);
    trendmap_rel = trendmap_thismonth ./ chlmean_thismonth .*100;
    trendmap_rel_stack(:,:,i) = trendmap_rel;
    pval_stack(:,:,i) = pval;
end

%% Plot
figure('position',[292    71   578   876]);
clf
for i = 1:length(month_ind)
    trendmap_rel = trendmap_rel_stack(:,:,i);
    pval = pval_stack(:,:,i);
    sigpts_m = pval < 0.05;
    haxes(i) = subplot(4,2,i);
    m_proj('lambert','long',[-80 -55],'lat',[-70 -60]); hold on
    m_pcolor(lon_chl,lat_chl,trendmap_rel); shading flat;
    colormap(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([-10 10]);
    [xs,ys] = m_ll2xy(lon_chl,lat_chl);
   stipple(xs,ys,sigpts_m,'markersize',6);
    if i == 8
        cb = colorbar;
        cb.Label.String = {'Relative Trend in Chl-a (% yr^{-1})';'::: p < 0.05'};
        cb.LineWidth = 1.5;
        cb.Location = 'southoutside';
        cb.Position = [0.275 0.08 0.5 0.025];
        cb.FontSize = 14;
        cb.Label.FontSize = 16;
    end

    % Add coast
    coastcolor = [.8 .8 .8];
    for ii=1:length(coastlon)
        p_coast = m_patch(coastlon{ii},coastlat{ii},coastcolor,'LineWidth',1);
    end; hold on;
    % Add ice shelf
    shelfcolor = [1 1 1];
    for ii=1:length(shlflon)
        p_shelf = m_patch(shlflon{ii},shlflat{ii},shelfcolor,'LineWidth',1);
    end
    % Add continental shelf break
    m_plot(contshlf_lon,contshlf_lat,'linestyle',':','color','k','linewidth',1.5);

    title(mostr{i},'fontsize',14,'fontweight','normal');
    hold off
    pos = get(gca,'position');
    set(gca,'position',[pos(1)-0.02,pos(2)+0.01,pos(3)+0.055,pos(4)+0.035]);
end

sublabel(haxes,1,14,'fontsize',16);

%% Save
print(gcf, 'Fig_9_Overall_Chl_Trends_ByMonth_Maps_Relative.png','-dpng','-r200');
print(gcf, 'Fig_9_Overall_Chl_Trends_ByMonth_Maps_Relative.eps','-depsc');

