PRECIPITATION ALGORITHM FOR MS04313

easypre.m

%%% EASYPRE.m

%%% EASYPRE (run in lowercase) is a script to take precipitation data and
%%% run it through a state algorithm. It is based on the differences
%%% between raw measurements. The output is a cell array called, usefully,
%%% "Results", which contains the following:

%%% DATENUM | ACCUMULATED RAIN | RAW GAUGE | BASELINE MEASUREMENTS |
%%% DIFFERENCE BETWEEN SUBSEQUENT MEASUREMENTS | FLAG

%%% See the state loop below for details.
%%% The initial measurement of the RawGauge is assumed to the the baseline,
%%% and it is likely not '0'. For accumulation beginning at 0, subtract
%%% this measurement universally from the column2 of accumulated
%%% measurement.

%%% This is currently running on the CENMET precipitation data from
%%% 2012-2014, and on the "SA" column. The column or reference can easily
%%% be changed to run on the "SH", or the two could be run in tandem.

%%% Output files can be generated by enabling the fopen commands.

%%% filename = 'output_CENMET_precipitation.csv';
%%% fid = fopen(filename,'w');
%%% fprintf(fid,'%s, %s, %s, %s, %s,
%%% %s\n','DATE','ACCUM','RAW','BASE','DIFF','FLAG');

%%% import the concatenated CENMET data. 
d = importdata('ALL_CEN_PRE.csv');

%%% Create a Results vector to hold ALL THE RECORDS! (Minus 1?)
Results = cell(length(d)-1,6);

%%% first two records
FirstLine = regexp(d{1},',','split');
SecondLine = regexp(d{2},',','split');

%%% Raw Gauge == Accumulated Gauge == Baseline
RawGauge = FirstLine{1,12}; AccGauge = FirstLine{1,12}; Baseline = FirstLine{1,12};
%%% RawGauge ?= Accumulated Gauge ?= Baseline
RawGauge2 = SecondLine{1,12}; AccGauge2 = SecondLine{1,12}; Baseline2 = SecondLine{1,12};

%%% Datestamps for records 1 and 2
DN = str2num(FirstLine{1,2});
DN2 = str2num(SecondLine{1,2});

%%% Turn off warning about converting strings to doubles versus to numbers
warning off

%%% The difference between the raw record and the accumulated record for
%%% the first Difference, always accepted
Diff = str2num(RawGauge) - str2num(AccGauge); Flag = 'A';
%%% The difference between subsequent raw records for the second record;
%%% will be tested
Diff2 = str2num(RawGauge2) - str2num(RawGauge);
FirstRecord = {DN, str2num(AccGauge), str2num(RawGauge), str2num(Baseline), Diff, Flag};

%%% States based on Diff

%%% If the difference is positive or nothing, the baseline rises,
%%% accumulation rises, and the data is accepted.
if Diff2 >= 0;
    Baseline2 = str2num(Baseline) + Diff2;
    AccGauge2 = str2num(AccGauge) + Diff2;
    Flag2 = 'A';
    
%%% If the difference is negative and it is much greater (1 order of
%%% magnitude give or take 0.005 to deal with the case of 0) that the
%%% previous difference, then the Baseline is reset to the raw gauge and it
%%% is flagged as a maintenance, no accumulation occurs. The flag is 'N'
%%% for 'mainteNance' because well, 'M' is already taken!


%%% NOTE: started with a factor of 0.005. If the diff is -0.1, though,
%%% then the log scaled diff was still 0.04 which was too large for
%%% that bound. 
elseif Diff2 < 0 & (log10(-1*Diff2+1) >= (abs(Diff) + 0.05));
    Baseline2 = str2num(RawGauge2);
    Flag2 = 'N';
    AccGauge2 = str2num(AccGauge) + 0;
    
%%% If the differences is negative and it is smaller than the previous
%%% difference after the aforementioned one-order-of-magnitude-ish buffer
%%% is applied, then the Baseline is adjusted to this negative difference,
%%% but the accumulation remains, and this baseline is used for the next
%%% measurement. The flag is 'Q' for questionable (Evaporation?)
elseif Diff2 < 0 & (log10(-1*Diff2+1) < abs(Diff) + 0.05);
    Baseline2 = str2num(Baseline) + Diff2;
    Flag2 = 'Q';
    AccGauge2 = str2num(AccGauge) + 0;
    
%%% If all these cases fail (for example a NaN in the data, which is the
%%% fail case I can think of), then the Flag is 'M' and the data is just
%%% repeated from the previous good measurement to keep the program
%%% running.
else Baseline2 = Baseline; Flag3 = 'M'; AccGauge2 = AccGauge; RawGauge2 = RawGauge;
end

%%% The second record reflects these criteria
SecondRecord = {DN2, AccGauge2, str2num(RawGauge2), Baseline2, Diff2, Flag2};

%%% Both the second and first record are given to the results.
Results(1,:) = FirstRecord;
Results(2,:) = SecondRecord;

%%%-- FOR FILE IO---%%%
%%% fprintf(fid,'%f, %f, %f, %f, %f, %s\n', DN, AccGauge, str2num(RawGauge), Baseline, Diff, Flag);
%%% fprintf(fid,'%f, %f, %f, %f, %f, %s\n', DN2, AccGauge2, str2num(RawGauge2), Baseline2, Diff2, Flag2);
%%%-- BACK TO REGULAR PROGRAM --- %%%

%%% A similar protocol to the above criteria is implemented in a loop, but
%%% in a cleaner, more efficient way which only passes in one row at a
%%% time.
for i = 3:length(Results);
    
    %%% The third line is parsed
    ThirdLine = regexp(d{i},',','split');
    
    %%% Raw measurement is taken from the third line
    RawGauge3  = ThirdLine{1,12};
    %%% Comparative measurements are taken from the Results array,
    %%% so they are all numeric (see else clause below!)
    CompareRaw = Results{i-1,3}; CompareAcc = Results{i-1,2};
    CompareBase = Results{i-1,4}; CompareDiff = Results{i-1,5};

    %%% The third DN comes from the raw data
    DN3 = str2num(ThirdLine{1,2});
    %%% The difference is the difference between the raw measurement 
    %%% and the running "baseline" which is not the same thing as the raw
    %%% or the accumulated (represents the 'state of the gauge')
    Diff3 = str2num(RawGauge3) - CompareBase;

    %%% similar to above, if the Diff > 0 or = 0 then the base and
    %%% accumulation reflect this change and it is accepted
    if Diff3 >= 0;
        Baseline3 = CompareBase + Diff3;
        Flag3 = 'A';
        AccGauge3 = CompareAcc + Diff3;
    
    %%% if the Diff is less than 0 by a whole lot, it is probably a
    %%% maintenance event, and the Flag is 'N' for maintenance, and the
    %%% Baseline is the raw, and the Accumulation does not accrue
    
    %%% NOTE: started with a factor of 0.005. If the diff is -0.1, though,
    %%% then the log scaled diff was still 0.04 which was too large for
    %%% that bound. 
    elseif Diff3 < 0 & (log10(-1*Diff3 + 1) > abs(CompareDiff) + 0.05);
        Baseline3 = CompareRaw;
        Flag3 = 'N';
        AccGauge3 = CompareAcc;
    
    %%% if the difference is less than 0 and the log of this difference is
    %%% less than the abs of the comparative difference plus a small factor
    %%% for the case of 0, then it's just regular old evaporation, and gets
    %%% a 'Q'
    elseif Diff3 < 0 & (log10(-1*Diff3 + 1) < abs(CompareDiff) + 0.05);
    
        Baseline3 = CompareBase + Diff3;
        Flag3 = 'Q';
        AccGauge3 = CompareAcc;
    
    %%% In the else clause, turn the CompareRaw into a string since
    %%% normally it would be numeric and the write out function reads in 
    %%% type string
    else Baseline3 = CompareBase; Flag3 = 'M'; AccGauge3 = CompareAcc; RawGauge3 = num2str(CompareRaw);
    end

    ThirdRecord = {DN3, AccGauge3, str2num(RawGauge3), Baseline3, Diff3, Flag3};    

    Results(i,:) = ThirdRecord;
    
    clearvars ThirdRecord
    %%%--- FOR FILE IO---%%%
    %%% fprintf(fid,'%f, %f, %f, %f, %f, %s\n',DN3, AccGauge3, str2num(RawGauge3), Baseline3, Diff3, Flag3)
    %%% --- BACK TO REGULAR PROGRAM---%%%
end

%%%--- FOR FILE IO---%%%
%%% fclose(fid);
%%%--- FOR FILE IO---%%%

fprintf(1,'%s\n','Hey, it's done');