path = '/Users/iljagolikov/Documents/ssh/study_Ilja/baka/period/fsum_icevol_'

pikkused = [];
aastad = 1993:1:2020;
algus= [];

lopp=[];
andmed = [];

for aasta = string(1993:1:2020)

ncf = append(path,aasta,'.nc');

data = ncread(ncf, 'sithick');
data = data(:);

st = find(data,1,'first');
ed = find(data,1,'last');

if sum(data)~=0
    data(1:st)= NaN;
    data(end:-1:ed)=NaN;
    pikkused(end+1) = length(data) - sum(isnan(data)); 
    algus(end+1)=st+1;
    lopp(end+1)=ed-1;
    andmed = cat(2,andmed,data); 
else
    pikkused(end+1) = NaN;     %kui polnud jaad siis pikkus = 0
    data(1:243)=0;           % andmed  = NaN
    algus(end+1)=NaN;         % algus = NaN
    lopp(end+1)=NaN;           % lopp = NaN
    andmed = cat(2,andmed,data);
end



k = plot(data)
hold on
k.LineWidth = 0.8
title('Kogu jää ruumala 1993/94-2020/21')
xlabel('Päevade arv alates 1-st oktoobrist') 
ylabel('Jää ruumala [km^3]') 
%grid on
disp(aasta);
end


%%
pikkused = pikkused(:);
aastad = aastad(:);
algus=algus(:);
lopp=lopp(:);

tabel = table(aastad,algus,lopp,pikkused);

%% Trendid

lm1 = fitlm(aastad,algus)
figure(1)
plot(lm1)

title('aastad/algus')

lm2 = fitlm(aastad,lopp)
figure(2)
plot(lm2)
title('aastad/lopp')

lm3 = fitlm(aastad,pikkused)
figure(3)
plot(lm3)
title('aastad/pikkused')

% Kogu jaapindala aastas
sum_ruumala = nansum(andmed)
lm4 = fitlm(aastad,sum_ruumala)
figure(4)
plot(lm4)
title('aastad/ice volume')

%% Statistika tabel
stats = ["min";"max";"mean";"std";"trend";"p-value"];

algus_stat = [min(algus),max(algus),nanmean(algus),nanstd(algus),table2array(lm1.Coefficients(2,1)),table2array(lm1.Coefficients(2,4))];
algus_stat=algus_stat(:);
lopp_stat = [min(lopp),max(lopp),nanmean(lopp),nanstd(lopp),table2array(lm2.Coefficients(2,1)),table2array(lm2.Coefficients(2,4))];
lopp_stat=lopp_stat(:);
pikkused_stat = [min(pikkused),max(pikkused),mean(pikkused),std(pikkused),table2array(lm3.Coefficients(2,1)),table2array(lm3.Coefficients(2,4))];
pikkused_stat=pikkused_stat(:);

% mean paevade kaupa ja std keskmistest paevade vaartustest
ice_vol = [min(min(andmed)),max(max(andmed)),nanmean(nanmean(andmed,2)),nanstd(nanmean(andmed,2)),table2array(lm4.Coefficients(2,1)),table2array(lm4.Coefficients(2,4))]
ice_vol=ice_vol(:)

stat = table(stats,algus_stat,lopp_stat,pikkused_stat,ice_vol)


%% main plot

% %keskmine jaa pindala paevas ule koikide aastate ja std-d
% mean_fsum = nanmean(andmed,2);
% std = nanstd(mean_fsum);
% std_plus = mean_fsum + std;
% std_minus = mean_fsum - std;
% 
% % min algus/lopp ja mean algus/lopp +- std plot
% plot(min(algus),0,'.','Color','k','MarkerSize', 20)
% hold on
% plot(max(lopp),0,'.','Color','k','MarkerSize', 20)
% hold on
% errorbar(nanmean(algus),0,nanstd(algus),'horizontal','Color','k','Marker','*','LineWidth',1,"CapSize",20)
% hold on
% errorbar(nanmean(lopp),0,nanstd(lopp),'horizontal','Color','k','Marker','*','LineWidth',1,"CapSize",20)
% hold on
% 
% %keskmise ja std plottimine
% plot(mean_fsum,'k','LineWidth',1)
% hold on
% x = 1:243;
% y1 = transpose(std_plus);
% y2 = transpose(std_minus);
% fill([x(49:219) fliplr(x(49:219))], [y1(49:219) fliplr(y2(49:219))],...
%     'b','FaceAlpha',0.2,'EdgeColor','k')
% 
% grid on; hold on
% 
% legend('start','end','mean start +-std',...
%     'mean end +- std','mean ice vol', 'mean ice vol +-std','Location','northwest')
% title('Kogu jää ruumala 1993/94-2020/21')
% xlabel('Päevade arv alates 1-st oktoobrist') 
% ylabel('Jää ruumala [km^3]') 
% %plot(data)
% axis([0 250 -1 2.5])




%% Koik ruumalad 2D
f = figure(1);
f.Position = [100 100 1200 500];
time_num = [];
for aasta = 1993:2020;
    for aeg = 1:243;
        time_num(end+1)=datenum(aasta,10,aeg) ;

    end
end

test = datestr(time_num());
time = datetime(test)
p = plot(time,andmed(:))
p.LineWidth = 1.5
%title('Kogu jää ruumala 1993/94-2020/21')

ax = gca;
ax.XAxis.FontSize = 16;
ax.YAxis.FontSize = 16;
ax.LineWidth = 1.2;

xlabel('Aastad','fontweight','bold','FontSize', 20) 
ylabel('Jää ruumala [km^3]','fontweight','bold','FontSize', 20)
%set(gca,'ytick',0:3000:15000);
datetick('x','yy')
grid on
%%

paevad = 1:243;
f = figure(1);
f.Position = [100 100 1000 800];

% miks ta ei joonesta 2019 ja 2020
fake = zeros(243,1)

%data = cat(2,data,fake);  % Kui runnid esimest korda votta komeaat maha
years = 1993:2021;

years = years(:);
data(data == 0) = NaN

time_num = [];
for aeg = 1:243;
        time_num(end+1)=datenum(1993,10,aeg) ;
end
test = datestr(time_num());
time = datetime(test)

pcolor(time,years,transpose(data));
shading flat
cb=colorbar;
grid on
cb.Label.String = 'Jää ruumala [km^3]'
cb.Label.FontWeight = 'bold'
cb.FontSize = 20
set(gca,'ytick',aastad(1:end))
%set(gca,'xtick',paevad(1:20:end))
%axis([0 243 1993 2021 0 15000])
datetick('x','mmm')
%title('Summaarne jää pindala 1993/94-2020/21')
ax = gca;
ax.XAxis.FontSize = 16;
ax.YAxis.FontSize = 16;
ax.LineWidth = 1;
xlabel('Kuud','fontweight','bold','FontSize', 20) 
ylabel('Aastad','fontweight','bold','FontSize', 20)









