% Explore CMEMS data - GlobColour
% Sources:
% https://resources.marine.copernicus.eu/product-detail/OCEANCOLOUR_GLO_CHL_L4_REP_OBSERVATIONS_009_082/DATA-ACCESS
% https://resources.marine.copernicus.eu/product-detail/OCEANCOLOUR_GLO_OPTICS_L4_REP_OBSERVATIONS_009_081/DATA-ACCESS
%
% Jessie Turner January 2022
% Updated to plot a time series of L4 chl from daily interpolated product
%
% Updated to compare daily to monthly
% By taking monthly average of all daily at each location and comparing to
% the given "monthly" file.

clear; close all hidden; format shortg;

% Set folder where data are located:
folder = '/Users/jsturner/satellite_data/WAP_CMEMS/';
% Get monthly file to use as comparison
file1 = [folder 'dataset-oc-glo-bio-multi-l4-chl_4km_monthly-rep_1641577754305.nc'];
chl_m0 = ncread(file1,'CHL');
monthly_time = ncread(file1,'time');
[~,monthly_date]=timesince1900(monthly_time);
chl_m = chl_4thorderpolynomial(chl_m0); % correct for factor of 2.

% load chlorophyll data coordinates lat, lon
load('/Users/jsturner/satellite_data/matlab_WAP/daily_CMEMS/Mean_Chl_20231026.mat')

% Load ALL Daily gap-filled chlorophyll (combined into one .nc file in
% previous version of this script)
folder2 = '/Users/jsturner/satellite_data/matlab_WAP/daily_CMEMS/';
load([folder2 'daily_chl_CMEMS_WAP_v4.mat']); % This is a LARGE file!
daily_date = date_mat; clearvars date_mat monthly_time;

%%
close all hidden;
% reshape into location x time 2D array
monthly_chl = reshape(chl_m,[(size(chl_m,1).*size(chl_m,2)),length(monthly_date)]);
daily_chl = reshape(chl,[(size(chl,1).*size(chl,2)),length(daily_date)]);

% for each month, average daily data from each location
monthly_mean_of_daily = nan(size(monthly_chl,1),size(monthly_chl,2));
for i = 1:length(monthly_date)
    thismonth = monthly_date(i);
    if i==length(monthly_date) % for last one, reference end of month manually
        nextmonth = datenum(2021,07,01);
    else
        nextmonth = monthly_date(i+1);
    end
    ind_daily = find(daily_date >= thismonth & daily_date < nextmonth);
    chl_thismonth = daily_chl(:, ind_daily);
    chl_tosave = nanmean(chl_thismonth,2);
    monthly_mean_of_daily(:,i) = chl_tosave;
end % ends month-to-month loop of 286 months total

% reshape into Nx1 arrays for Scatter Plot colored by Kernel Density
% Estimate (KDE) by Nils "scatter_kde.m"
mm_chl = reshape(monthly_chl,[numel(monthly_chl),1]);
md_chl = reshape(monthly_mean_of_daily,[numel(monthly_chl),1]);

clearvars chl_thismonth chl_tosave daily_chl

% Plot
disp(' ');
disp('Calculating Probability density estimate of points for figure...');
figure
clf
% plot(mm_chl,md_chl,'ko');
% scatter_kde(mm_chl,md_chl,'MarkerSize',100);
binscatter(mm_chl,md_chl,250);
hold on
xx = 0.01:0.01:40; yy=xx; plot(xx,yy,'k-','linewidth',1.5); % black line 1:1
makepretty;
xlabel('Monthly Chl-a (mg m^{-3})');
ylabel('Monthly Mean of Daily Gap-Filled Chl-a (mg m^{-3})');
set(gca,'yscale','log','xscale','log');
cb = colorbar;
cb.Label.String = 'Probability Density Estimate';
cb.LineWidth = 1.5;
cb.Label.FontSize = 14;
caxis([0 10000]); 
colormap(turbo);
axis([0.05 50 0.05 50]);
hold off

% Calculate general statistics
% uses function "quickstats.m"
% First, remove nans
md_chl(isnan(mm_chl))=[]; mm_chl(isnan(mm_chl))=[];
mm_chl(isnan(md_chl))=[]; md_chl(isnan(md_chl))=[];
[bias,mean_ratio, MAE, RMSE,meanAPD,R] = quickstats(md_chl,mm_chl)

annotation('textbox',...
    [0.15 0.65 0.25 0.25],...
    'String',{['Bias = ' num2str(round(bias,2))],...
    ['Mean Ratio = ' num2str(round(mean_ratio,2))],...
    ['MAE = ' num2str(round(MAE,2))],...
    ['RMSE = ' num2str(round(RMSE,2))],...
    ['meanAPD = ' num2str(round(meanAPD,2))],...
    ['R = ' num2str(round(R,2))]},...
    'FontSize',14,'EdgeColor','k','LineWidth',2,...
    'BackgroundColor','w','Color','k');


% Save figure
print(gcf,'Fig_S3_Compare_GapFilled_NonGapFilled.png','-dpng','-r300');
print(gcf,'Fig_S3_Compare_GapFilled_NonGapFilled.eps','-depsc');



