%dataset use example
%author: Oleksandr (Sasha) Huziy


of1 = importdata('test_1yr/MONTREAL.OF1');
headers = cell(length(of1.colheaders),1);
i = 1;
%delete stars from column names
for name = of1.colheaders
   tempo = name{1};
   tempo(strfind(tempo, '*')) = '';
   headers{i} = tempo;
   i = i + 1;
end

%starting date of your data
start_date = datenum('1960/01/01');

%% create dataset
columns = cell(length(headers), 1);

for col = 1: length(columns)
    columns{col,1} = of1.data(:,col);
end

ds = dataset(columns{1:end}, 'VarNames', headers);

%create date column
ds.DATES = arrayfun(@(d) datevec(addtodate(start_date, d, 'day')), ds.DAY - 1, 'UniformOutput', false);
%create month column
ds.MONTH = arrayfun(@(the_date) the_date{1}(2) , ds.DATES);
ds.MONTH_DAY = arrayfun(@(the_date) the_date{1}(3), ds.DATES);

stamp_year = 1961;


%calculate group stats daily means
stats = grpstats(ds, {'MONTH','MONTH_DAY'}, {'mean', 'min', 'max'}, 'DataVars', {'QE','QH'});
stats.STAMP_DATES = arrayfun(@(m,d) datenum(stamp_year, m, d), stats.MONTH, stats.MONTH_DAY);


summary(stats);
%plot results
figure
    hold on
    plot(stats.STAMP_DATES, stats.mean_QE, 'r');
    plot(stats.STAMP_DATES, stats.mean_QH, 'g');
    
    title('Daily climatology')
    ylabel('W/m**2')
    legend('QE', 'QH');
    
    grid on
    
    %put months on x-axis
    datetick('x','dd-mmm','keepticks')
    hold off


figure
    hold on
    

    
    X = [ stats.STAMP_DATES', fliplr(stats.STAMP_DATES')];
    Y = [ stats.max_QE', fliplr(stats.min_QE')];
    
       
    hf = fill(X, Y, 'g');
    
    h1 = plot(stats.STAMP_DATES, stats.min_QE, 'b');
    h2 = plot(stats.STAMP_DATES, stats.max_QE, 'r');

    
    %set xaxis ticks and limits
    tick_vals = arrayfun(@(month) datenum(stamp_year, month, 15), 1:2:12);
    set(gca,'XTick',tick_vals)
    
    %put months on x-axis
    datetick('x','dd-mmm', 'keepticks');
   
    xlim([stats.STAMP_DATES(1) stats.STAMP_DATES(end)])
   
    
    title('Ranges of daily means')
    ylabel('W/m**2')
    
    
    
    
    legend([h1, h2], 'QE(min)', 'QE(max)');
    legend boxoff
    legend boxon
    grid on
    
    hold off


%display(stats);

%annual means
stats = grpstats(ds, {'YEAR'}, {'mean', 'max', 'min'}, 'DataVars', {'QE','QH', 'SWE'});

%linear tendency
p = polyfit(stats.YEAR, stats.mean_QE, 1);
line_data_vals = p(1) * stats.YEAR + p(2);
figure
    hold on
    
    
    
    title('Annual means')
    ylabel('W/m**2')
    h1 = plot(stats.YEAR, stats.max_QE, 'k--');
    h2 = plot(stats.YEAR, stats.mean_QE, 'b');
    h3 = plot(stats.YEAR, line_data_vals, 'b--');
    plot(stats.YEAR, stats.min_QE,'k--')
    
    
    legend([h1 h2 h3],'min/max', 'mean', 'linear tend')
    grid on
    
    hold off

%linear tendency
p = polyfit(stats.YEAR, stats.mean_QE, 1);

%%Plot swe
f = figure
    hold on
    
    h = plot(ds.SWE)
    grid on
    print(f, 'swe.png', '-dpng')

    hold off





