% Plot each month mean of Chl-a  over all years.
% May 17, 2022 this was Figure S2.
% Figure 3 as of August 17, 2022.
% Figure 2 as of October 23, 2023.

%% Load data
datafolder = '/Users/jsturner/satellite_data/matlab_WAP/daily_CMEMS/';
load([datafolder 'daily_chl_CMEMS_WAP_v3.mat']);
% Load monthly chl to use dates to help generate monthly averages:
folder = '/Users/jsturner/satellite_data/WAP_CMEMS/';
file1 = [folder 'dataset-oc-glo-bio-multi-l4-chl_4km_monthly-rep_1641577754305.nc'];
chl_m = ncread(file1,'CHL');
monthly_time = ncread(file1,'time');
[~,monthly_date]=timesince1900(monthly_time);
% Load lat lon from monthly chl
lon0 = ncread(file1,'lon');
lat0 = ncread(file1,'lat');
lon_chl = repmat(lon0,1,size(chl,2));
lat_chl = repmat(lat0,1,size(chl,1))';
clearvars *0;
% 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:
coastcolor = [1 1 1];
load([mapfolder 'continental_shelf_WAP.mat']);
% Set map axes limits
axislims=[-79.9792  -55.0208  -69.9792  -60.0208];
load('Calgae2.mat');

%% Trim daily chl dataset spatially to not include ice-covered locations
% Calculate how many points are present for each pixel
time_len = size(chl,3); % number of timesteps (>8000 days)
perc_cov = sum(isfinite(chl),3)./time_len;
% Trim to only analyze pixels that have >30% of the record.
chl_data = reshape(chl,[(size(perc_cov,1).*size(perc_cov,2)),time_len]);
chl_data((perc_cov <= 0.3),:) = nan;
perc_cov_data = reshape(perc_cov,(size(perc_cov,1).*size(perc_cov,2)),1);
perc_cov_data(perc_cov <= 0.3) = 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)]);

%% Organize data into monthly means
daily_date = date_mat; clearvars date_mat monthly_time;
monthly_chl = reshape(chl_m,[(size(chl_m,1).*size(chl_m,2)),length(monthly_date)]);
daily_chl = reshape(chl_trimmed,[(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
md_chl = reshape(monthly_mean_of_daily,[size(chl,1),size(chl,2),size(chl_m,3)]);
clearvars folder chl_thismonth chl_tosave daily_chl daily_date ans chl_m daily_date err file* i ind_daily monthly_chl monthly_mean_of_daily nextmonth thismonth

%% Get mean and median of each month over all years
yrlist = 1998:2022; % think of it being the year it ends on (Feb)
date_columns = datevec(monthly_date);
yyyy = date_columns(:,1);
mm = date_columns(:,2);
% Pre-allocate
monthclim_chl = [];
monthclim_med_chl = [];
for k = 1:12
    ind_thismo = find(mm==k);
    chl_thismo = squeeze(md_chl(:,:,ind_thismo));
    % how many points are included (data availability?)
    % Calculate how many points are present for each pixel
    time_len = size(chl_thismo,3); % number of timesteps
    perc_cov = sum(isfinite(chl_thismo),3)./time_len;
    % Trim to only analyze pixels that have >60% of the years for that month.
    chl_data = reshape(chl_thismo,[(size(perc_cov,1).*size(perc_cov,2)),time_len]);
    chl_data((perc_cov <= 0.75),:) = nan;
    perc_cov_data = reshape(perc_cov,(size(perc_cov,1).*size(perc_cov,2)),1);
    perc_cov_data(perc_cov <= 0.75) = nan;
    chl_trimmed = reshape(chl_data,[size(perc_cov,1),size(perc_cov,2),time_len]);
    
    mean_thismo = nanmean(chl_trimmed,3);
    median_thismo = nanmedian(chl_trimmed,3);
    monthclim_chl = cat(3,monthclim_chl,mean_thismo);
    monthclim_med_chl = cat(3,monthclim_med_chl,median_thismo);
end
clearvars *_thismo k mm

%% Plot
close all hidden;
mostr = {'July',...
    'August','September','October','November','December',...
    'January','February','March','April','May','June',};

figure('position',[292   108   678   900]);
clf
mo_ind = [7:12,1:6];
for k = 1:12
    mm = mo_ind(k);
    haxes(k) = subplot(4,3,k);
    m_proj('lambert','long',[-80 -55],'lat',[-70 -60]); hold on
    m_pcolor(lon_chl,lat_chl,squeeze(monthclim_chl(:,:,mm)));
    m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,'fontsize',14,'linewidth',1,'fontname','Arial','backcolor','none','xticklabels',[],'yticklabels',[]);
    shading flat;
    caxis([0.2 4]);
    colormap(Calgae2)
    set(gca,'colorscale','log','ytick',[],'xtick',[]);
    
    if k==11
        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,11); 
        cw.TickLabels = num2cell(round(logspace(-2,1,11),1));
        cw.LineWidth = 1.5;
        cw.Position = [0.25   0.08   0.55    0.025];
        cw.Label.FontSize = 14;
        cw.FontSize = 14;
    end
    title(char(mostr(k)),'fontweight','normal','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',1.5);
    % scooch subplots slightly
    pos = get(gca,'position');
    set(gca,'position',[pos(1)-0.02,pos(2)+0.01,pos(3)+0.04,pos(4)+0.04]);
end

sublabel(haxes,12,8,'fontsize',14);
print(gcf,'Fig_2_Month_Clim_All.png','-dpng','-r200');
print(gcf,'Fig_2_Month_Clim_All.eps','-depsc');
