% Plot sea ice date of advance and retreat 1997-2020
% Using data from Sharon Stammerjohn, CUBoulder
%
% Info from Sharon:
%
% Note, these image files are gridded into ~25km by ~25km grid cells and
% includes the entire Southern Ocean (with continent and islands masked
% out).  The grid cell areas are approximate because the grid is in
% polarstereographic projection and is only equal-area at latitude 70S,
% with grid cells becoming smaller/bigger north/south of 70S.  If you need
% to know the actual area of the grid cells, I can include that info too,
% as well as the lat/lon info.
% 
% I didn?t include all of the geospatial ?coordinate reference systems?
% (crs) info in that *.nc file, but you can find that info in this file,
% which is just an example of another product using this same exact
% projection and crs info:
% "NSIDC0081_SEAICE_PS_S25km_20220711_v2.0.nc"
%  
% Code by Jessie Turner July 14, 2022
% Adapted December 2023 using correct Lat, Lon
clear; close all hidden;

% Load map information
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 'continental_shelf_WAP.mat']);

% Load sea ice dates data from Sharon
folder = '/Users/jsturner/satellite_data/WAP_Ice_Stammerjohn/';
file1 = 'iceret_9720.nc';
file2 = 'iceadv_9720.nc';

lat = ncread([folder file1], 'Longitude'); % actually latitude
lon = ncread([folder file1], 'Latitude'); % actually longitude in units of degrees 0 to 360
lon(lon>=180) = lon(lon>=180)-360; % get longitude in correct units

lat_wap = lat(35:97, 145:235);
lon_wap = lon(35:97, 145:235);

%% Load data and organize
% Load: Date of sea ice retreat %%%%%%%%%%%%%%%%
% Load all years and plot average
yearlist = (1997:2020)'; % year day ranging up to 410 (mid Feb of year 2)
yrstrlist = num2str(yearlist);
ice_ret = nan(size(lat,1),size(lat,2),length(yearlist));
for i = 1:length(yearlist)
    eval(['str1 = "iceret_' yrstrlist(i,:) '";']);
    thisyr = ncread([folder file1],str1);
    ice_ret(:,:,i) = thisyr;
%     figure;pcolor(thisyr); shading flat; caxis([180 500]);
end
% Fill non-data pixels with NaN
ice_ret(ice_ret==0) = NaN;
ice_ret(ice_ret==-999) = NaN;
% Calculate finite number of measurements and get rid of points without
% sufficient data over the time series:
time_len = size(ice_ret,3); % number of timesteps (>8000 days)
perc_cov = sum(isfinite(ice_ret),3)./time_len;
% % Trim to only analyze pixels that have >80% of the record.
acceptable = .8;
ret_data = reshape(ice_ret,[(size(perc_cov,1).*size(perc_cov,2)),time_len]);
ret_data((perc_cov <= acceptable),:) = nan;
perc_cov_data = reshape(perc_cov,(size(perc_cov,1).*size(perc_cov,2)),1);
perc_cov_data(perc_cov <= acceptable) = nan;
ret_trimmed = reshape(ret_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)]);
perc_cov_wap = perc_cov(35:97, 145:235);
clearvars icedata perc_cov

% Take mean
mean_ice_ret = nanmean(ret_trimmed,3);
mean_ice_ret_wap = mean_ice_ret(35:97, 145:235);

% Load: Date of sea ice advance %%%%%%%%%%%%%%
% Load all years and plot average
ice_adv = nan(size(lat,1),size(lat,2),length(yearlist));
for i = 1:length(yearlist)
    eval(['str1 = "iceadv_' yrstrlist(i,:) '";']);
    thisyr = ncread([folder file2],str1);
    ice_adv(:,:,i) = thisyr;
%     figure;pcolor(thisyr); shading flat; caxis([180 500]);
end
% Fill non-data pixels with NaN
ice_adv(ice_adv==0) = NaN;
ice_adv(ice_adv==-999) = NaN;
% Calculate finite number of measurements and get rid of points without
% sufficient data over the time series:
time_len = size(ice_adv,3); % number of timesteps (>8000 days)
perc_cov = sum(isfinite(ice_adv),3)./time_len;
% % Trim to only analyze pixels that have >80% of the record.
adv_data = reshape(ice_adv,[(size(perc_cov,1).*size(perc_cov,2)),time_len]);
adv_data((perc_cov <= acceptable),:) = nan;
perc_cov_data = reshape(perc_cov,(size(perc_cov,1).*size(perc_cov,2)),1);
perc_cov_data(perc_cov <= acceptable) = nan;
adv_trimmed = reshape(adv_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)]);
perc_cov_wap = perc_cov(35:97, 145:235);

% Take mean
mean_ice_adv = nanmean(adv_trimmed,3);
mean_ice_adv_wap = mean_ice_adv(35:97, 145:235);


%% Calculate trends
% Calculate trend in retreat
t1 = (1:length(yearlist))';
disp('Calculating trend at each x,y point (may be slow)...');
trendmap_raw = nan(length(ret_data),1);
pval_raw = nan(length(ret_data),1);
for ii = 1:size(ret_data,1)
    data = [t1,ret_data(ii,:)'];
    [m,b] = TheilSen(data);
    line0 = t1.* m + b;
    % [~,pval0] = corr(data); % test whether the time axis and the data are correlated
    [~,pval0] = corr([line0 ret_data(ii,:)']); % test whether the line and the data are correlated
    trendmap_raw(ii) = m;
    pval_raw(ii) = pval0(2);
end
trendmap_ret = reshape(trendmap_raw,[size(perc_cov,1),size(perc_cov,2)]);
pval_ret = reshape(pval_raw,[size(perc_cov,1),size(perc_cov,2)]);
% [trendmap_ret,pval_ret] = trend(ret_trimmed,t1,'omitnan','Type','Kendall');
% [~,pval_ret] = mann_kendall(ret_trimmed);
sigpts_ret = pval_ret < 0.05;
trendmap_ret_wap = trendmap_ret(35:97, 145:235);
pval_ret_wap = pval_ret(35:97, 145:235);
sigpts_ret_wap = sigpts_ret(35:97, 145:235);

% Calculate trend in ice advance date
trendmap_raw = nan(length(adv_data),1);
pval_raw = nan(length(adv_data),1);
for ii = 1:size(adv_data,1)
    data = [t1,adv_data(ii,:)'];
    [m,b] = TheilSen(data);
    line0 = t1.* m + b;
    % [~,pval0] = corr(data); % test whether the time axis and the data are correlated
    [~,pval0] = corr([line0 adv_data(ii,:)']); % test whether the line and the data are correlated
    trendmap_raw(ii) = m;
    pval_raw(ii) = pval0(2);
end
trendmap_adv = reshape(trendmap_raw,[size(perc_cov,1),size(perc_cov,2)]);
pval_adv = reshape(pval_raw,[size(perc_cov,1),size(perc_cov,2)]);
% [trendmap_adv,p] = trend(adv_trimmed,t1,'omitnan','Type','Kendall');
% [~,pval_adv] = mann_kendall(adv_trimmed);
sigpts_adv = pval_adv < 0.05;
trendmap_adv_wap = trendmap_adv(35:97, 145:235);
pval_adv_wap = pval_adv(35:97, 145:235);
sigpts_adv_wap = sigpts_adv(35:97, 145:235);


%% Figures: Date of sea ice retreat
close all hidden;
% Plot amount of data
figure('position',[100 100 850 730]);
clf;

% Plot mean ice retreat date
plt(1) = subplot(2,2,1);
m_proj('lambert','long',[-80 -55],'lat',[-70 -60]); hold on
m_pcolor(lon_wap,lat_wap,mean_ice_ret_wap);shading flat;
m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,'fontsize',14,'linewidth',1,'fontname','Arial','backcolor','none','xticklabels',[],'yticklabels',[]);
colormap(plt(1),turbo(25));
caxis([244 425]);
cb = colorbar;
dates_ticks = [213 244 274 305 335 366 397 425 456];
cb.Ticks = dates_ticks;
% cb.TickLabels = list_doy_strings(dates_ticks);
cb.TickLabels = {'S','O','N','D','J','F','M'};
cb.FontSize = 14;
cb.LineWidth = 1.5;
cb.Location = 'southoutside';
cb.Label.String = 'Mean Day-of-Retreat';
cb.Label.FontSize = 16;
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; hold on;
% 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);
pos = get(gca,'position');
set(gca,'pos',[pos(1)-0.01 pos(2)-0.01 pos(3)+0.05 pos(4)+0.05]);
    

% Plot mean ice advance date
plt(2) = subplot(2,2,2);
m_proj('lambert','long',[-80 -55],'lat',[-70 -60]); hold on
m_pcolor(lon_wap,lat_wap,mean_ice_adv_wap);shading flat;
m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,'fontsize',14,'linewidth',1,'fontname','Arial','backcolor','none','xticklabels',[],'yticklabels',[]);
colormap(plt(2),turbo(25));
caxis([32 244]);
cb = colorbar;
% dates_ticks = [213 244 274 305 335 366 397 425 456];
dates_ticks = [32 60 91 121 152 182 213 244];
cb.Ticks = dates_ticks;
% cb.TickLabels = list_doy_strings(dates_ticks);
cb.TickLabels = {'F','M','A','M','J','J','A','S'};
cb.FontSize = 14;
cb.LineWidth = 1.5;
cb.Location = 'southoutside';
cb.Label.String = 'Mean Day-of-Advance';
cb.Label.FontSize = 16;
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; hold on;
% 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);
pos = get(gca,'position');
set(gca,'pos',[pos(1)-0.01 pos(2)-0.01 pos(3)+0.05 pos(4)+0.05]);
    
% Plot trend in ice retreat date
plt(3)=subplot(2,2,3);
m_proj('lambert','long',[-80 -55],'lat',[-70 -60]); hold on
m_pcolor(lon_wap,lat_wap,trendmap_ret_wap); shading flat;
m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,'fontsize',14,'linewidth',1,'fontname','Arial','backcolor','none','xticklabels',[],'yticklabels',[]);
colormap(plt(3),flipud(cbrewer('div','RdBu',31)));
caxis([-4 4]);
hold on
% [xs,ys] = m_ll2xy(lon_wap,lat_wap);
cb = colorbar;
cb.Label.String = {'Trend in day of retreat (days yr^{-1})';'::: p < 0.05'};
cb.LineWidth = 1;
cb.Location = 'southoutside';
cb.FontSize = 14;
cb.Label.FontSize = 16;
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; hold on;
% 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);
pos = get(gca,'position');
set(gca,'pos',[pos(1)-0.01 pos(2)-0.01 pos(3)+0.05 pos(4)+0.05]);
m_plot(lon_wap(sigpts_ret_wap),lat_wap(sigpts_ret_wap),'k.','markersize',12);

% Plot trend in ice advance date
plt(4) = subplot(2,2,4);
m_proj('lambert','long',[-80 -55],'lat',[-70 -60]); hold on
m_pcolor(lon_wap,lat_wap,trendmap_adv_wap); shading flat;
m_grid('tickdir','in','xtick',-80:10:-55,'ytick',-70:2:-60,'fontsize',14,'linewidth',1,'fontname','Arial','backcolor','none','xticklabels',[],'yticklabels',[]);
colormap(plt(4),flipud(cbrewer('div','RdBu',31)));
caxis([-4 4]);
hold on
% [xs,ys] = m_ll2xy(lon_wap,lat_wap);
% plot(xs,ys,sigpts_adv_wap,'k.','markersize',18);
cb = colorbar;
cb.Label.String = {'Trend in day of advance (days yr^{-1})';'::: p < 0.05'};
cb.LineWidth = 1;
cb.Location = 'southoutside';
cb.FontSize = 14;
cb.Label.FontSize = 16;
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; hold on;
% 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);
pos = get(gca,'position');
set(gca,'pos',[pos(1)-0.01 pos(2)-0.01 pos(3)+0.05 pos(4)+0.05]);
m_plot(lon_wap(sigpts_adv_wap),lat_wap(sigpts_adv_wap),'k.','markersize',12);


sublabel(plt,1,32,'fontsize',18);
% % save figure
print(gcf,'Fig_S16_SeaIce_RetreatAdvance.png','-dpng','-r300');

%%
% Display information about statistics
pval_ret = pval_ret_wap(lon_wap > -80 & lon_wap < -55 & lat_wap > -70 & lat_wap < -60);
pval_adv = pval_adv_wap(lon_wap > -80 & lon_wap < -55 & lat_wap > -70 & lat_wap < -60);

pval_ret_sig = pval_ret(pval_ret<0.05);
pval_adv_sig = pval_adv(pval_adv< 0.05);
pval_ret_notsig = pval_ret(pval_ret > 0.05);
pval_adv_notsig = pval_adv(pval_adv > 0.05);

disp(['Ice Retreat Date Trend significant at ' ...
    num2str(numel(pval_ret_sig)) ' x,y locations, p < ' num2str(nanmax(pval_ret_sig))]);
disp(['Ice Retreat Date Trend not sig. at ' ...
    num2str(2081 - numel(pval_ret_sig)) ' x,y locations, p > ' num2str(nanmin(pval_ret_notsig))]);
disp(['Ice Advance Date Trend significant at ' ...
    num2str(numel(pval_adv_sig)) ' x,y locations, p < ' num2str(nanmax(pval_adv_sig))]);
disp(['Ice Advance Date Trend not sig. at ' ...
    num2str(2081 - numel(pval_adv_sig)) ' x,y locations, p > ' num2str(nanmin(pval_adv_notsig))]);

%%
% hold on;
% for i=1:length(coastlon)
%     p_coast = patch(coastlon{i},coastlat{i},coastcolor,'LineWidth',1);
% end
% % Add ice shelf
% shelfcolor = [1 1 1];
% for i=1:length(shlflon)
%     p_shelf = patch(shlflon{i},shlflat{i},shelfcolor,'LineWidth',1);
% end
