% Compile daily record of SST for 1997 to 2022 from CMEMS SST
% Data Source:
% Product identifier
% SST_GLO_SST_L4_REP_OBSERVATIONS_010_011
% Product name
% Global Ocean OSTIA Sea Surface Temperature and Sea Ice Reprocessed
% Dataset
% Single dataset in this product:
% METOFFICE-GLO-SST-L4-REP-OBS-SST
%
% https://data.marine.copernicus.eu/product/SST_GLO_SST_L4_REP_OBSERVATIONS_010_011/services
%
% Downloaded in February 2023
% Downloaded data from June 1, 1997 to May 31, 2022 in 7 files.
%
% Merged into one time series using script
% "calculate_daily_SST_timeseries_WAP.m"

clear; close all hidden; format shortg;

% Load data
datafolder = '/Users/jsturner/satellite_data/matlab_WAP/SST_data_Mar2023/';
load([datafolder 'daily_SST_WAP.mat']);
disp(' ');
disp('File loaded.');

% Convert from Kelvin to Celsius
sst = sst - 273.15;

% Load map information AND REGIONS
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']);

% Trim to only analyze pixels that have >95% of the record.
time_len = size(sst,3); % number of timesteps
perc_cov = sum(isfinite(sst),3)./time_len;
perc_perc = perc_cov.*100;
sst_data = reshape(sst,[(size(perc_cov,1).*size(perc_cov,2)),time_len]);
sst_data((perc_cov <= 0.95),:) = nan;
perc_cov_data = reshape(perc_cov,(size(perc_cov,1).*size(perc_cov,2)),1);
perc_cov_data(perc_cov <= 0.95) = nan;
sst_trimmed = reshape(sst_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)]);

%% Snip to each decade and establish polygons for regions

% Decade 1 = July 1, 2000 to June 30, 2010 (00-01 to 09-10)
% Decade 2 = July 1, 2010 to June 30, 2020 (10-11 to 19-20)
% ind_0010 = find(date_mat >= datenum(2000,07,01) & date_mat <= datenum(2010,06,30));
% ind_1020 = find(date_mat >= datenum(2010,07,01) & date_mat <= datenum(2020,06,30));
ind_0010 = find(date_mat >= datenum(2000,07,01) & date_mat <= datenum(2011,06,30));
ind_1020 = find(date_mat >= datenum(2011,07,01) & date_mat <= datenum(2022,06,30));
date_0010 = date_mat(ind_0010);
date_1020 = date_mat(ind_1020);
sst_0010 = sst_trimmed(:,:,ind_0010);
sst_1020 = sst_trimmed(:,:,ind_1020);
nyrs = 11;

% Set up spatial indexes within the lon,lat of the sst data
spf = inpolygon(lon_sst,lat_sst,lon_spf,lat_spf);
miz = inpolygon(lon_sst,lat_sst,lon_miz,lat_miz);
nos = inpolygon(lon_sst,lat_sst,lon_nos,lat_nos);
mis = inpolygon(lon_sst,lat_sst,lon_mis,lat_mis);
sos = inpolygon(lon_sst,lat_sst,lon_sos,lat_sos);

disp('Split into 2 decades, indices established for region polygons.'); 

%% Daily - Loop by region
close all hidden;

% For each region
regionlist = {'spf','miz','nos','mis','sos'};
titlelist = {'Southern Polar Front','Marginal Ice Zone','Northern Shelf',...
    'Mid Shelf','Southern Shelf'};

% Start figure
figure('position',[100   100   846   816]);
clf
ylimlist = [-1 8; -3 4; -3 4; -3 4; -3 4];

for i = 1:length(regionlist)
    % i = 1;
    reg1 = char(regionlist(i));
    disp(char(titlelist(i)));
    % calculate number of pixels in this polygon
    eval(['no_pts = length(find(' reg1 '==1));']);

    % preallocate and fill a per-day sst plot of whole time series
    sst_0010_reg = nan(size(sst_0010,3),1);
    sst_1020_reg = nan(size(sst_1020,3),1);
    % spatial avg per each day of time series
    for ii = 1:size(sst_0010,3)
        sstthisday = sst_0010(:,:,ii);
        eval(['sstthisday_reg = nanmean(sstthisday(' reg1 '));']);
        sst_0010_reg(ii) =  sstthisday_reg;
    end
    for ii = 1:size(sst_1020,3)
        sstthisday = sst_1020(:,:,ii);
        eval(['sstthisday_reg = nanmean(sstthisday(' reg1 '));']);
        sst_1020_reg(ii) =  sstthisday_reg;
    end

    % reshape by day of year
    doy_0010 = nan(length(date_0010),1);
    doy_1020 = nan(length(date_1020),1);
    for ii = 1:length(date_0010)
        doy_0010(ii) = date2jd(date_0010(ii)) - 0.5;
    end
    for ii = 1:length(date_1020)
        doy_1020(ii) = date2jd(date_1020(ii)) - 0.5;
    end
    sst_0010_reg_stack = nan(nyrs,366); % 10 years by 366 days
    sst_1020_reg_stack = nan(nyrs,366);
    for ii = 1:366
        ind_thisdoy = find(doy_0010==ii);    
        sst_thisdoy = sst_0010_reg(ind_thisdoy);
        if size(sst_thisdoy,1)==nyrs
            sst_0010_reg_stack(:,ii) = sst_thisdoy;
        elseif size(sst_thisdoy,1)==nyrs-1
            sst_0010_reg_stack(:,ii) = [sst_thisdoy;nan];
        end
    end
    for ii = 1:366
        ind_thisdoy = find(doy_1020==ii);    
        sst_thisdoy = sst_1020_reg(ind_thisdoy);
        if size(sst_thisdoy,1)==nyrs
            sst_1020_reg_stack(:,ii) = sst_thisdoy;
        elseif size(sst_thisdoy,1)==nyrs-1
            sst_1020_reg_stack(:,ii) = [sst_thisdoy;nan];
        end
    end

    % change days <180 to be >366 for continuity
    seas_0010 = [sst_0010_reg_stack(:,181:366),sst_0010_reg_stack(:,1:180)];
    seas_1020 = [sst_1020_reg_stack(:,181:366),sst_1020_reg_stack(:,1:180)];
    x = 181:546;

    % Statistics
    % test whether decades are significantly different
    pval = nan(366,1);
    for ii = 1:366
        testdata = [seas_0010(:,ii),seas_1020(:,ii)];
        if any(isnan(testdata(:)))
            pval(ii) = NaN;
        elseif round(mean(testdata(:,1))) == round(mean(testdata(:,2)))
            pval(ii) = NaN;
        else
        [~,~,stats] = kruskalwallis(testdata,[1 2],'off');
        [c,~,~,~] = multcompare(stats,'estimate','kruskalwallis','display','off');
        pval(ii) = c(6);
        end
    end
    sig = find(pval < 0.05);
    x_sig = x(sig);
    disp(['Min p value: ' num2str(nanmin(pval))]);
    disp(['Max p value: ' num2str(nanmax(pval))]);
    
    % take mean, median, stdev, and percentiles of each decade
    mean_seas_0010 = nanmean(seas_0010,1);
    mean_seas_1020 = nanmean(seas_1020,1);
    median_seas_0010 = nanmedian(seas_0010,1);
    median_seas_1020 = nanmedian(seas_1020,1);
    
    % Plots
    lb = [198,219,239]./255; % light blue
    
    haxes(i) = subplot(3,2,i);
    plot(x,seas_0010,'.','color',lb);
    hold on
    plot(x,seas_1020,'.','color',lb);
    p1 = plot(x,mean_seas_0010,'b-','linewidth',2);
    p2 = plot(x,mean_seas_1020,'m-','linewidth',2);
    title(titlelist(i))
    makepretty;
    ylim(ylimlist(i,:));
    dates_ticks = [213 244 274 305 335 366 397 425 456 486];
    datestr_bi = list_doy_strings(dates_ticks);
    xticks(dates_ticks);
    set(gca,'xminortick','off');
%     xticklabels(datestr_bi);
    xticklabels({'Aug','Sep','Oct','Nov','Dec','Jan','Feb','Mar','Apr','May'});
    haxes(i).XTickLabelRotation = 90;
    xlim([195 510])
% %     plot x locations in season that decades are sig different
%     yc = ylimlist(i,1) + 0.75;
%     plot(x_sig,ones(length(x_sig),1).*yc,'k.','markersize',14);

    if i==1  
        ylabel('SST (^o C)');
    elseif i==3
        ylabel('SST (^o C)');
    elseif i==5
        ylabel('SST (^o C)');
    end

    if i==1
        ll = legend([p1,p2],'2001-2011','2012-2022');
        ll.Location = 'northeast';
        ll.FontSize = 16;
    end
     
end
sublabel(haxes,15,20,'fontsize',18);

%%%%%%%%%%%% Add map as 6th subplot %%%%%%%%%%%%%%%%%%
mean_sst = nanmean(sst_trimmed,3);
subplot(3,2,6);
cla;
m_proj('lambert','long',[-80 -55],'lat',[-70 -60]); hold on
m_pcolor(lon_sst,lat_sst,mean_sst); shading flat;
m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,...
    'fontsize',14,'linewidth',1,'fontname','Arial','backcolor',...
    'none','xticklabels',[],'yticklabels',[]);
caxis([-2 5]);
colormap(turbo(35))
cw = colorbar;
cw.LineWidth = 1;
cw.Label.String = 'Mean SST (^o C)';
cw.Location = 'southoutside';
cw.FontSize = 14;
cw.Label.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
% Plot regions %%% NEW as of July 15, 2022
plt1 = m_plot(lon_sos,lat_sos,'-','linewidth',2);
plt2 = m_plot(lon_mis,lat_mis,'-','linewidth',2);
plt3 = m_plot(lon_nos,lat_nos,'-','linewidth',2);
plt4 = m_plot(lon_miz,lat_miz,'-','linewidth',2);
plt5 =  m_plot(lon_spf,lat_spf,'-','linewidth',2);
str1 = {'a','b','c','d','e'};
str_x = [-77.69,-77.69,-61.4,-70.23,-74.9];
str_y = [-61.6,-66,-63.23,-66.35,-68.56];
m_text(str_x,str_y,str1,'fontsize',18, 'fontweight','bold');
% Realign and Label
pos = get(gca,'position');
set(gca,'position',[pos(1)-0.04,pos(2)-0.06,pos(3)+0.06,pos(4)+0.06]);


% Save Figure
print(gcf, 'Fig_S14_SST_Decade_Seasonal_Regions.png','-dpng','-r200');


