% 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 May 2022
% Goal: plot trends from Monthly data including maps and region timeseries
% Version 2 August 2023, to address reviewer comments using Mann-Kendall.
% Version 3 October 2023, de-seasonalized before calculating trends.

clear; close all hidden;

%% Load data
% Load monthly chl file 
datafolder = '/Users/jsturner/satellite_data/matlab_WAP/means_CMEMS/';
load([datafolder 'Month_Means_v3.mat']);
datafolder2 = '/Users/jsturner/satellite_data/matlab_WAP/daily_CMEMS/';
load([datafolder2 'Mean_Chl_20231026.mat']);
load('Calgae2.mat');

%% 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 'WAP_regions_2022-07-15.mat']);
load([mapfolder 'continental_shelf_WAP.mat']);

%% Clean up data
 % Remove pixels that have < 70% 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.7),:) = nan;
perc_cov_data = reshape(perc_cov,(size(perc_cov,1).*size(perc_cov,2)),1);
perc_cov_data(perc_cov <= 0.7) = nan;
chl_trimmed = reshape(chl_data,[size(perc_cov,1),size(perc_cov,2),time_len]);
perc_cov_trimmed = reshape(perc_cov_data,[size(perc_cov,1),size(perc_cov,2)]);

%% Average per region, plot time series, save statistics
close all hidden
reg_str_list = {'spf','miz','nos','mis','sos'};
title_list = {'Southern Polar Front','Marginal Ice Zone','Northern Shelf','Mid Shelf','Southern Shelf'};
date_monthly = (datenum(date_columns));
% Set up figure
figure('position',[377   442   863   536])
clf

% Preallocate arrays for slope, R, and p-value
trends_reg = nan(length(reg_str_list),1);
trends_reg_rel = nan(length(reg_str_list),1);
p_val_reg = nan(length(reg_str_list),1);
for j = 1:length(reg_str_list)
    disp(' ');
    disp(title_list(j));
    % Index region
    reg_str = char(reg_str_list(j));
    eval(['ind_reg = inpolygon(lon_chl,lat_chl,lon_' reg_str ',lat_' reg_str ');']);

    % Preallocate arrays
    chl_reg = nan(size(md_chl,3),1);
    
     % Average values for each month and save to array
    for i = 1:size(md_chl,3)
        chl_thismo = squeeze(md_chl(:,:,i));
        chl_reg(i) = nanmean(chl_thismo(ind_reg)); 
    end
    
    % Long-term mean chl for this region
    chl_reg_avg = nanmean(chl_reg);
    
    % Deal with NaN - clip them out
    date_clipped = date_monthly; date_clipped(isnan(chl_reg))=[];
    chl_clipped = chl_reg; chl_clipped(isnan(chl_reg))=[];
%     % Deal with NaN - fill via interpolation
%     date_clipped = date_monthly;
%     chl_clipped = fillmissing(chl_reg,'linear','SamplePoints',date_monthly);  

    % Detrend by seasonal cycle ("deseason" Climate Data Toolbox)
    % This function assumes
    %    y = y_0 + y_tr + y_season + y_var + y_noise
    % where
    % y_0 is the long-term mean,
    % y_trend is the long-term trend,
    % y_season is the typical seasonal anomaly, which the season function obtains after detrending and removing the mean,
    % y_var represents interannual variability, and
    % y_noise is everything else
    %
    % The deseason function simply works like this:
    % Ads = A - season(A,t);
    chl_reg_ds = deseason(chl_clipped,date_clipped);
    
    % Plot time series at that location
    haxes(j) = subplot(2,3,j);
    makepretty;
    hold on
    plot(date_monthly,chl_reg,'.-','linewidth',1.5,'markersize',12);% ,'markersize',6);
    plot(date_clipped,chl_reg_ds,'.-','linewidth',1.5,'markersize',12);
    title(title_list(j));
    datetick('x','keeplimits');
    xlim([date_monthly(1)-120 date_monthly(end)+120]);
    if j==1
        ylabel('Chl-a (mg m^{-3})');
        ll = legend('Chl-a','De-seasonalized');
    elseif j==4
        ylabel('Chl-a (mg m^{-3})');
    end
    ylim([0.08 5]);
    set(gca,'yscale','log','ytick',[0.1 0.25 1 2.5 5],'yticklabels',{'0.1','0.25' '1' '2.5' '5'});
    
%     % Calculate trend
%     [m,b,R2,eq1,p_val,f] = lineofbestfit(date_clipped,chl_reg_ds,1,0);
% 
%     %%%%%%%%%% NEW: CALCULATE MANN-KENDALL TREND
%     datain = [date_clipped,chl_reg_ds];alpha = 0.05; wantplot = 0;
%     [taub,tau,~,sig,Z,~,~,~,~,~,~,~,~,~,~]= ktaub(datain, alpha, wantplot);
%     disp(['p = ' num2str(sig)]);
% 
%     % Save to arrays
%     trends_reg(j) = m;
%     p_val_reg(j) = sig;
%     tr = m.*12; % Per month to Per year
%     rel_tr = tr ./ chl_reg_avg .* 100; % Percent, relative to long-term mean chl for that region.
%     trends_reg_rel(j) = rel_tr; % Percent, relative to long-term mean chl for that region.
%     disp(['Trend = ' num2str(tr) ' mg m-3 yr-1']);
%     disp(['Relative trend = ' num2str(rel_tr) ' % yr-1']);
end

%%%%%%%%%%%% Add map as 6th subplot %%%%%%%%%%%%%%%%%%
subplot(2,3,6);
cla;
m_proj('lambert','long',[-80 -55],'lat',[-70 -60]); hold on
m_pcolor(lon_chl,lat_chl,chl_mean); shading flat;
m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,...
    'fontsize',14,'linewidth',2,'fontname','Arial','backcolor',...
    'none','xticklabels',[],'yticklabels',[]);
caxis([0.2 4]);
colormap(Calgae2)
set(gca,'colorscale','log');
cw = colorbar;
cw.LineWidth = 1;
cw.Label.String = 'Chl-a (mg m^{-3})';
cw.Location = 'southoutside';
cw.Ruler.Scale = 'log'; c1.Ruler.MinorTick = 'on';
cw.Ticks = logspace(-2,1,13); 
cw.TickLabels = num2cell(round(logspace(-2,1,13),1));
cw.FontSize = 14;
cw.Label.FontSize = 14;
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
% 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',2);

% Plot regions as of July 15, 2022
plt1 = m_plot(lon_sos,lat_sos,'-','linewidth',2);
plt2 = m_plot(lon_mis,lat_mis,'-','linewidth',2);
plt3 = m_plot(lon_nos,lat_nos,'-','linewidth',2);
plt4 = m_plot(lon_miz,lat_miz,'-','linewidth',2);
plt5 =  m_plot(lon_spf,lat_spf,'-','linewidth',2);
str1 = {'a','b','c','d','e'};
str_x = [-77.69,-77.69,-61.4,-70.23,-74.9];
str_y = [-61.6,-66,-63.23,-66.35,-68.56];
m_text(str_x,str_y,str1,'fontsize',18, 'fontweight','bold','color','w');
% Realign
pos = get(gca,'position');
set(gca,'position',[pos(1)-0.03,pos(2)-0.03,pos(3)+0.06,pos(4)+0.06]);

sublabel(haxes);

% % save figure
set(gcf, 'InvertHardCopy', 'off');
print(gcf,'Fig_S10_Chl_Trends_ByRegion.png','-dpng','-r200');
% 
