% 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.
% To create only Table 4, October 2023, for MEPS submission.
% Version 3 October 2023, de-seasonlizes Chl-a 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']);

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

%% 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);
Sp_rho_reg = 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))=[];

    % 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);
    
    % Calculate Trend:
    %%%%%%%%%% Linear Regression
%     [m,b,R2,eq1,p_val,f] = lineofbestfit(date_clipped,chl_reg_ds,1,0);
    %%%%%%%%%% Theil-Sen Slope non-parametric
    data = [date_clipped,chl_reg_ds];
    [m,b] = TheilSen(data);
%     line0 = date_clipped .* m + b;
%     [~,pval0] = corr([line0,chl_reg_ds]);
    %%%%%%%%%% Mann-Kendall non-parametric
%     datain = [date_clipped,chl_reg_ds];alpha = 0.05; wantplot = 0;
%     [taub,tau,~,sig,Z,~,~,~,~,~,~,~,~,~,~]= ktaub(datain, alpha, wantplot);
%     disp(['p = ' num2str(sig)]);
    [~,pval]=mann_kendall(chl_reg_ds);
    disp(['p = ' num2str(pval)]);

    % Save to arrays
%     p_val_reg(j) = pval0(2);
    p_val_reg(j) = pval;
    tr = m.*12; % Per month to Per year
    trends_reg(j) = tr;
    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

%% Create table
close all hidden;
format shortg;

trends_str = cellstr(num2str(round(trends_reg,6),'%.5f'));
t1 = table(title_list',trends_str,round(trends_reg_rel,3),round(p_val_reg,4));
var_names_t1 = {'Region','Trend_peryr','Relative_Trend_peryr','p_value'};
t1.Properties.VariableNames = var_names_t1;
writetable(t1,'Table_4_overall_Chl_trends_byRegion_20240209.csv','WriteVariableNames',true);


