%%%% FIGURE 3 for Postdoc Paper 1
% Explore CMEMS data - GlobColour
% Sources:
% https://resources.marine.copernicus.eu/product-detail/OCEANCOLOUR_GLO_CHL_L4_REP_OBSERVATIONS_009_082/DATA-ACCESS
% As of September 7, 2022

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 avg chl
load([datafolder 'Mean_Chl_20231026.mat']);
% Load map information
mapfolder = '/Users/jsturner/satellite_data/matlab_WAP/WAP_map_SchuylerJack/';
load([mapfolder 'WAP_regions_2022-07-15.mat']);


%% Average over all years for all bloom indices
% Long-term means of all indices
bi_avg_thres = nanmean(bi_thres,3);

% Fix Tmax to be 365 days and more rather than 1-180.
tmax(tmax<180) = tmax(tmax<180)+365;
tmax_avg = nanmean(tmax,3);

max_chl_avg = nanmean(max_chl,3);

%% Plot each bloom index long-term average start date (map)
close all hidden;

%% Calculate average indices for each region
% Index regions
spf = inpolygon(lon_chl,lat_chl,lon_spf,lat_spf);
miz = inpolygon(lon_chl,lat_chl,lon_miz,lat_miz);
nos = inpolygon(lon_chl,lat_chl,lon_nos,lat_nos);
mis = inpolygon(lon_chl,lat_chl,lon_mis,lat_mis);
sos = inpolygon(lon_chl,lat_chl,lon_sos,lat_sos);

% Start date
meanStart_spf = nanmean(bi_avg_thres(spf));
meanStart_miz = nanmean(bi_avg_thres(miz));
meanStart_nos = nanmean(bi_avg_thres(nos));
meanStart_mis = nanmean(bi_avg_thres(mis));
meanStart_sos = nanmean(bi_avg_thres(sos));
errStart_spf = nanstd(bi_avg_thres(spf));
errStart_miz = nanstd(bi_avg_thres(miz));
errStart_nos = nanstd(bi_avg_thres(nos));
errStart_mis = nanstd(bi_avg_thres(mis));
errStart_sos = nanstd(bi_avg_thres(sos));

startdate_jd = round([meanStart_spf;meanStart_miz;meanStart_nos;meanStart_mis;meanStart_sos]);
startdate_str = list_doy_strings(startdate_jd);
startdate_err_days = round([errStart_spf;errStart_miz;errStart_nos;errStart_mis;errStart_sos]);

% Tmax
meanTmax_spf = nanmean(tmax_avg(spf));
meanTmax_miz = nanmean(tmax_avg(miz));
meanTmax_nos = nanmean(tmax_avg(nos));
meanTmax_mis = nanmean(tmax_avg(mis));
meanTmax_sos = nanmean(tmax_avg(sos));
errTmax_spf = nanstd(tmax_avg(spf));
errTmax_miz = nanstd(tmax_avg(miz));
errTmax_nos = nanstd(tmax_avg(nos));
errTmax_mis = nanstd(tmax_avg(mis));
errTmax_sos = nanstd(tmax_avg(sos));

Tmax_jd = round([meanTmax_spf;meanTmax_miz;meanTmax_nos;meanTmax_mis;meanTmax_sos]);
Tmax_str = list_doy_strings(Tmax_jd);
Tmax_err_days = round([errTmax_spf;errTmax_miz;errTmax_nos;errTmax_mis;errTmax_sos]);

%%%%%%%%%% Mean Chl
meanchl_spf = nanmean(chl_mean(spf));
meanchl_miz = nanmean(chl_mean(miz));
meanchl_nos = nanmean(chl_mean(nos));
meanchl_mis = nanmean(chl_mean(mis));
meanchl_sos = nanmean(chl_mean(sos));
errchl_spf = nanstd(chl_mean(spf));
errchl_miz = nanstd(chl_mean(miz));
errchl_nos = nanstd(chl_mean(nos));
errchl_mis = nanstd(chl_mean(mis));
errchl_sos = nanstd(chl_mean(sos));

MeanChl = round([meanchl_spf;meanchl_miz;meanchl_nos;meanchl_mis;meanchl_sos],2);
MeanChlErr = round([errchl_spf;errchl_miz;errchl_nos;errchl_mis;errchl_sos],2);

% Max Chl
meanMaxC_spf = nanmean(max_chl_avg(spf));
meanMaxC_miz = nanmean(max_chl_avg(miz));
meanMaxC_nos = nanmean(max_chl_avg(nos));
meanMaxC_mis = nanmean(max_chl_avg(mis));
meanMaxC_sos = nanmean(max_chl_avg(sos));
errMaxC_spf = nanstd(max_chl_avg(spf));
errMaxC_miz = nanstd(max_chl_avg(miz));
errMaxC_nos = nanstd(max_chl_avg(nos));
errMaxC_mis = nanstd(max_chl_avg(mis));
errMaxC_sos = nanstd(max_chl_avg(sos));

MaxChl = round([meanMaxC_spf;meanMaxC_miz;meanMaxC_nos;meanMaxC_mis;meanMaxC_sos],2);
MaxChl_err = round([errMaxC_spf;errMaxC_miz;errMaxC_nos;errMaxC_mis;errMaxC_sos],2);

%% Table
region_names = {'Southern_Polar_Front';'Marginal_Ice_Zone';'Northern_Shelf';...
    'Mid_Shelf';'Southern_Shelf'};
t1 = table(region_names,startdate_str,startdate_err_days,...
    Tmax_str,Tmax_err_days,MeanChl,MeanChlErr,MaxChl,MaxChl_err);
writetable(t1,'Table_1_Region_Avg_Bloom_Indices_20231027.csv','WriteVariableNames',true);



