eventCharacteristics.m 6.36 KB
Newer Older
Claas Faber's avatar
Claas Faber committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184
function [lengths,growths,decays,event_dates,peak_strengths,mean_strengths] = eventLengths(time_series,date_vector,threshold,type)

% Takes a time series of anomalies, date vector, and threshold and finds
% all events and their lengths in the time series for either positive or
% negative events. Event lengths are given in units of time step
% of the time axis, assuming that time stepping is regular.
%
% USAGE: [lengths,growths,decays,event_dates,peak_strengths,mean_strengths] = ...
%            eventLengths(time_series,date_vector,threshold,type)
%
% lengths, growths,
% decays, event_dates:      vectors that contain one value per event.
% time_series:              the data, a time series of anomalies
% threshold:                threshold above the mean (0 for anomalies) for which
%                           an anomaly is potentially part of an event
% type:                     1 for positive events, -1 for negative events
%
% An event is identified when anomalies of the same sign exceed a
% threshold with respect to the standard deviation of the anomaly type. E.g.,
% for positive events, take calculate the mean of all positive events
% and apply the threshold for identification to this mean,
% instead of to the standard deviation of the entire anomaly
% time series. (--> This is what eventLengths() does.)
%
% In comparison to eventLengths(), relativeEventLengths() sets
% the threshold for events differently and allows for
% systematic differences between positive and negative
% events. This is especially important when tackling
% issues of ENSO asymmetry.

 % ------------- %
 % QUALITY CHECK %
 % ------------- %

     % date vector and data of same size?
     n_ts   = length(time_series);
     n_time = size(date_vector,1);

     if ( n_ts ~= n_time )
         error('Dimension missmatch between time series and date vector.');
     end

% ------------ %
% PREPARE DATA %
% ------------ %

    % flip data upside down for negative events
    anomalies = type.*time_series;

    % mean of anomalies of interest, i.e. of either positive or
    % negative anomalies; for type = -1, the time series has
    % been flipped already
    % identify all desired anomalies
    is_anomaly = ( anomalies > 0 );
    relative_base      = std(anomalies(is_anomaly));

    % for verification purposes
    complementary_base = std(anomalies(~is_anomaly));

    % move offset according to threshold
    offset = anomalies-(threshold*relative_base);

    % date numbers for easier handling
    date_numbers = datenum(date_vector);

% --------------- %
% IDENTIFY EVENTS %
% --------------- %

    % set the time series to 1 for all positive events
    base_line = sign(offset);
    base_line(base_line<0) = 0;

    % breaks between events
    event_boundaries = diff(base_line);
    is_boundary = ( event_boundaries ~= 0 );

    % find indivudal events
    event_starts = find(event_boundaries > 0);
    event_stops  = find(event_boundaries < 0);

    % if there is either no event or no event_stops:
    % return NaN values
    if ( numel(event_starts) == 0 | numel(event_stops) == 0 )
        lengths     = NaN;
        growths     = NaN;
        decays      = NaN;
        event_dates = NaN;
        peak_strengths = NaN;
        mean_strengths = NaN;
        return;
    end

    % ignore partial events at the beginning nd end of the time series
    % leading events: ignore the first partial event
    if event_stops(1) < event_starts(1)
        event_stops = event_stops(2:end);
    end

    % trailing event: ignore the last (started) partial event
    if event_starts(end) > event_stops(end)
        event_starts = event_starts(1:end-1);
    end

    % add +1 to all event starts to make sure that they only
    % start when the anomaly value is already higher than the threshold
    event_starts = event_starts+1;

    % find events that last only for 1 month and ignore them
    micro_events = find( (event_stops-event_starts) < 2 );
    event_starts(micro_events) = [];
    event_stops(micro_events)  = [];

    % a check: there must be as many starting events as ending events
    if ( length(event_starts) ~= length(event_stops) )
        error('Missed leading or trailing events.');
    end

% -------------------------------- %
% IDENTIFY PEAKS AND COUNT LENGTHS %
% -------------------------------- %

    % number of events
    n_events = length(event_starts);

    % safe data
    lengths        = zeros(1,n_events);
    growths        = zeros(1,n_events);
    decays         = zeros(1,n_events);
    event_dates    = zeros(1,n_events);    % as datenums
    peak_strengths = zeros(1,n_events);
    mean_strengths = zeros(1,n_events);

    % loop over all events
    for c_event = 1:n_events

        % IDENTIFY EVENTS AND DIAGNOSE THEIR LENGTHS
        % ++++++++++++++++++++++++++++++++++++++++++

            % event cropping
            c_start    = event_starts(c_event);
            c_stop     = event_stops(c_event);
            event_data = anomalies(c_start:c_stop);

            % total event length
            event_length     = length(c_start:c_stop);
            lengths(c_event) = event_length;

            % growing and decaying phases (this is a bit tricky,
            % tread cautiously here)
            peak_position = find(event_data == max(event_data));

            % if more than one peak exists: good chance of corrupted data
            % ignore
            if ( length(peak_position) > 1 )
                lengths        = NaN;
                growths        = NaN;
                decays         = NaN;
                event_dates    = NaN;
                peak_strengths = NaN;
                mean_strengths = NaN;
                return;
            end

            % absolute position for the date vector
            abs_peak = c_start+peak_position-1;

            % save date of peak
            event_dates(c_event) = date_numbers(abs_peak);

            % growth and decay lengths in months
            % count peak-position twice, since growth ends there
            % and decay starts
            growths(c_event) = length(1:peak_position);
            decays(c_event)  = length(peak_position:event_length);

            % event strengths: peak strength and mean strength
            peak_strengths(c_event) = type*anomalies(abs_peak);   % to restore the correct sign
            mean_strengths(c_event) = type*mean(event_data);

    end

    % convert event peak date numbers into a more usable date vector
    event_dates = datevec(event_dates);