% Plot time series from regions of the WAP
% for NASA IDS WAP Project
% Jessie Turner
% February 2022
% Using CMEMS GlobColour Data, chl, daily-interpolated (level-4).
% Updated August 2022.
% Figures 9 through 12 of the paper as of August 20, 2022.
% Only Table 3 now in MEPS submission October 2023.

clear; close all hidden; format shortg;

years = 1998:2022;
austral_yr_str = {'97-98','98-99','99-00','00-01','01-02','02-03','03-04',...
        '04-05','05-06','06-07','07-08','08-09','09-10','10-11','11-12',...
        '12-13','13-14','14-15','15-16','16-17','17-18','18-19','19-20',...
        '20-21','21-22'};

% Load region polygons, map information, and latitude/longitude.
mapfolder = '/Users/jsturner/satellite_data/matlab_WAP/WAP_map_SchuylerJack/';
load([mapfolder 'WAP_regions_2022-07-15.mat']);
load([mapfolder 'bathy_and_coast.mat'],'coastlon','coastlat','shlflon','shlflat');
load('chldatacoords.mat');

% Load annual values of day-of-year of bloom initiation and other indices.
startdatefolder = '/Users/jsturner/satellite_data/matlab_WAP/daily_CMEMS/';
load([startdatefolder 'StartDates_CMEMS_WAP_v4.mat']);

% Change dates of tmax to be 180-400s in Julian day
tmax(tmax<180) = tmax(tmax<180)+365;

% Calculate days between bi and tmax for each location, each year.
bi = bi_thres;
tinc = tmax - bi; % in days.
tinc(tinc<0) = nan;

% Set up strings for region data and subplot titles for regions
reg_str_list = {'spf','miz','nos','mis','sos'};
title_list = {'Southern Polar Front','Marginal Ice Zone','Northern Shelf',...
    'Mid Shelf','Southern Shelf'};

ind_0010 = find(years >= 2001 & years <= 2011); % Jan 2001 to Jan 2010
ind_1020 = find(years >= 2011 & years <= 2021); % Jan 2011 to Jan 2020

%% For Each Region: Start Date
close all hidden

% Preallocate arrays for slope, R, and p-value
trends_reg = nan(length(reg_str_list),1);
R_val_reg = nan(length(reg_str_list),1);
p_val_reg = nan(length(reg_str_list),1);

for j = 1:length(reg_str_list)
    % 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
    bi_reg = nan(length(years),1);
     % Average values for each year and save to arrays
    for i = 1:length(years)
        bi_thisyr = squeeze(bi(:,:,i));
        bi_reg(i) = nanmean(bi_thisyr(ind_reg)); 
    end
    bi_reg_0010 = round(nanmedian(bi_reg(ind_0010)));
    bi_reg_1020 = round(nanmedian(bi_reg(ind_1020)));

    % Calculate trend
    date_clipped = years'; date_clipped(isnan(bi_reg))=[];
    bi_clipped = bi_reg; bi_clipped(isnan(bi_reg))=[];   
    %%%%%%%%%% Linear Regression  
%     [m,b,R,eq1,pval,f] = lineofbestfit(date_clipped,bi_clipped,1,0);
    %%%%%%%%%% Theil-Sen Slope non-parametric for trend
    data = [date_clipped,bi_clipped];
    [m,b] = TheilSen(data);
%     line0 = date_clipped .* m + b;
%     [~,pval0] = corr([line0,bi_clipped]);
    %%%%%%%%%% Mann-Kendall non-parametric for significance
    [~,pval]=mann_kendall(bi_clipped);
     % Save to arrays
    trends_reg(j) = round(m,3);
%     p_val_reg(j) = round(pval0(2),3);
    p_val_reg(j) = round(pval,3);
end
bi_trends_reg = trends_reg;
bi_pval_reg = p_val_reg;


%% For each region:  Peak Date
close all hidden;
clearvars trends_reg;

% Preallocate arrays for slope, R, and p-value
trends_reg = nan(length(reg_str_list),1);
p_val_reg = nan(length(reg_str_list),1);

for j = 1:length(reg_str_list)
    % 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
    tmax_reg = nan(length(years),1);
     % Average values for each year and save to arrays
    for i = 1:length(years)
        tmax_thisyr = squeeze(tmax(:,:,i));
        tmax_reg(i) = nanmean(tmax_thisyr(ind_reg)); 
    end
     % Calculate trend
    date_clipped = years'; date_clipped(isnan(tmax_reg))=[];
    tmax_clipped = tmax_reg; tmax_clipped(isnan(tmax_reg))=[];
        %%%%%%%%%% Linear Regression  
%     [m,b,R,eq1,pval,f] = lineofbestfit(date_clipped,bi_clipped,1,0);
    %%%%%%%%%% Theil-Sen Slope non-parametric for trend
    data = [date_clipped,tmax_clipped];
    [m,b] = TheilSen(data);
%     line0 = date_clipped .* m + b;
%     [~,pval0] = corr([line0,bi_clipped]);
    %%%%%%%%%% Mann-Kendall non-parametric for significance
    [~,pval]=mann_kendall(tmax_clipped);
     % Save to arrays
    trends_reg(j) = round(m,3);
%     p_val_reg(j) = round(pval0(2),3);
    p_val_reg(j) = round(pval,3);
end
tmax_trends_reg = trends_reg;
tmax_pval_reg = p_val_reg;

%% Clean data and save to spreadsheet
close all hidden;
t1 = table(title_list',bi_trends_reg,bi_pval_reg,...
    tmax_trends_reg,tmax_pval_reg);
t1.Properties.VariableNames = {'Ecoregion','Trend_StartDate_PerYear',...
    'p_value_StartDate','Trend_PeakDate_PerYear','p_value_PeakDate'};
writetable(t1,'Table_3_Trends_BloomIndices_20240209.csv','WriteVariableNames',true);


