% Convert waterbase data for vlugtenburg analyse
% Run oetsettings.m first

% check if you have ran oetsettings.m
if exist ('setProperty.m') ~= 2
    error ('DUDE...  run OETsettings first man!')
end


filename='id1-HOEKVHLD-200904010000-201112312359.txt'
% filename='id1-SCHEVNGN-200904010000-201012312359.txt'  %%% MISSCHIEN IS SCHEVNINGEN HANDIGER... MINDER RIVIER INPUT
% this file can be downloaded from waterbase via the following link
% http://live.waterbase.nl/wboutput.cfm?loc=HOEKVHLD&wbwns=1%7CWaterhoogte+in+cm+t.o.v.+normaal+amsterdams+peil+in+oppervlaktewater&byear=2009&bmonth=04&bday=01&eyear=2011&emonth=12&eday=31&output=Tekst&whichform=2

% Convert txt file to structure
DAT = rws_waterbase_read_edit(filename)

% save in handy matrix
h=[DAT.data.datenum DAT.data.waarde];
h_text='waterlevel @ HvH. (column 1 = Matlab time, column 2 = h in cm)';
cd ..
save ('HvHh.mat','h','h_text')

disp('HvHh.mat created')

%plot wave timeseries
figure
plot(DAT.data.datenum,DAT.data.waarde);
datetick('x')
title({ DAT.data.location ; DAT.data.waarnemingssoort})

%% count peaks
% first determine zero crossings.
% it is too difficutl for matlab to do this on the whole dataset so devide
% the dataset in subsets based on the length
peaks = [];

subsets = 100;

for k= 1:subsets
    if k==1
        sub_begin=1
        sub_end = floor(length(DAT.data.datenum)/subsets*k)
    else
        sub_begin=floor(length(DAT.data.datenum)/subsets*(k-1)+1)
        sub_end=floor(length(DAT.data.datenum)/subsets*k)
    end
    % determine zero crossings
    [xi,yi] = polyxpoly(DAT.data.datenum(sub_begin:sub_end),DAT.data.waarde(sub_begin:sub_end), ...
        DAT.data.datenum(sub_begin:sub_end),zeros(length((sub_begin:sub_end)),1));
    % extract the maximum values between the zero crossings
    for i = 1:(length(xi)-1)
        if abs(xi(i)-xi(i+1))>0.01
            raw_peaks(i) = max(DAT.data.waarde(DAT.data.datenum>xi(i) &DAT.data.datenum<xi(i+1)));
        end
    end

    peaks = [peaks raw_peaks(raw_peaks>1)];
    clear raw_peaks
end





