%% Kaart
file = 'ice_climatology5_2000_2020.nc';

lon = ncread(file,'lon');
lat = ncread(file,'lat');
avg_ice_days = ncread(file,'average_ice_days_per_year');

avg_ice_days(avg_ice_days < 1) = -1;   

figure('Color','w')

imagesc(lon,lat,avg_ice_days')
set(gca,'YDir','normal')
axis equal tight
box on

caxis([0 160])         
cm = parula(256);
cm = [0 0 0; cm];      
colormap(cm)
colorbar

xlabel('Longitude')
ylabel('Latitude')
title('Average ice days per year (2000–2020)')


%% Aegread
file = 'ice_climatology5_2000_2020.nc';

yearday = ncread(file,'yearday_range');

GoR_mean = ncread(file,'ice_extent_GoR_seasonal_mean');
GoR_std  = ncread(file,'ice_extent_GoR_seasonal_std');

GoF_mean = ncread(file,'ice_extent_GoF_seasonal_mean');
GoF_std  = ncread(file,'ice_extent_GoF_seasonal_std');

GoR_mean(isnan(GoR_mean)) = 0;  GoR_std(isnan(GoR_std)) = 0;
GoF_mean(isnan(GoF_mean)) = 0;  GoF_std(isnan(GoF_std)) = 0;

shift_day = 244; 

idx = yearday >= shift_day;

yearday_s = [yearday(idx)-365; yearday(~idx)];
GoR_mean_s = [GoR_mean(idx); GoR_mean(~idx)];
GoR_std_s  = [GoR_std(idx);  GoR_std(~idx)];

GoF_mean_s = [GoF_mean(idx); GoF_mean(~idx)];
GoF_std_s  = [GoF_std(idx);  GoF_std(~idx)];


figure('Color','w')

% GoR
subplot(2,1,1)
hold on
fill([yearday_s; flipud(yearday_s)], ...
     [GoR_mean_s-GoR_std_s; flipud(GoR_mean_s+GoR_std_s)], ...
     'b','FaceAlpha',0.25,'EdgeColor','none')
plot(yearday_s, GoR_mean_s,'b','LineWidth',2)
hold off
xticks([-122 -91 -61 -30 0 31 59 90 120 151 181 212])
xticklabels({'Sep','Oct','Nov','Dec','Jan','Feb','Mar','Apr','May','Jun','Jul','Aug'})
xlabel('Month')

grid on
ylim([-5000 max(GoR_mean_s+GoR_std_s)*1.05])
xlabel('Day of year (shifted)')
ylabel('Ice extent')
title('Gulf of Riga ice extent climatology')

% GoF
subplot(2,1,2)
hold on
fill([yearday_s; flipud(yearday_s)], ...
     [GoF_mean_s-GoF_std_s; flipud(GoF_mean_s+GoF_std_s)], ...
     'r','FaceAlpha',0.25,'EdgeColor','none')
plot(yearday_s, GoF_mean_s,'r','LineWidth',2)
hold off
xticks([-122 -91 -61 -30 0 31 59 90 120 151 181 212])
xticklabels({'Sep','Oct','Nov','Dec','Jan','Feb','Mar','Apr','May','Jun','Jul','Aug'})
xlabel('Month')

grid on
ylim([-1e4 max(GoF_mean_s+GoF_std_s)*1.05])
xlabel('Day of year (shifted)')
ylabel('Ice extent')
title('Gulf of Finland ice extent climatology')
