calcAnomalies.m 3.76 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
function [anomalies,cropped_dates] = calcAnomalies(data,date_vector,varargin)

% Calculates anomalies relative to a user-specified trend
% and to the climatological seasonal cycle.
%
% USAGE: [anomalies,cropped_dates] = calcAnomalies(data,date_vector)
%        [anomalies,cropped_dates] = calcAnomalies(__,trend)
%        [anomalies,cropped_dates] = calcAnomalies(__,start_year,stop_year);
%        [anomalies,cropped_dates] = calcAnomalies(__,trend,start_year,stop_year);
%
%        data and anomalies are vectors. date_vector is a Matlab date vector.
%        cropped_dates is the cropped date vector.
%
% Optional parameter are
%
%        trend:      trend that is to be subtracted from the data
%                    prior to the removal of the seasonal cycle.
%                    Vector that has the same length as data.
%                    Default: When no trend is provided, the linear
%                    trend is calculated and subtracted from
%                    the data.
%
%        start_year: When should anomaly calculation begin?
%
%        stop_year:  When should it stop?
%
% IF YOU DON'T WHISH TO CALCULATE A TREND,
% specify a "dummy trend" that consists of zeros and
% has the same size as the data.

% ----------- %
% PARSE INPUT %
% ----------- %

    n_additional = length(varargin);

    % define defaults:
        trend_data = [];
        start_year = min(date_vector(:,1));
        stop_year  = max(date_vector(:,1));

    switch n_additional
        case 0
        case 1
            trend_data = varargin{1};
        case 2
            start_year = varargin{1};
            stop_year  = varargin{2};
        case 3
            trend_data = varargin{1};
            start_year = varargin{2};
            stop_year  = varargin{3};
        otherwise
            error('Wrong amount of input variables specified.');
    end

% -------------- %
% CONSTRAIN TIME %
% -------------- %

    is_time =  ( date_vector(:,1) >= start_year ) & ( date_vector(:,1) <= stop_year );

    % cropping
    date_vector = date_vector(is_time,:);
    data        = data(is_time);

    % need to crop the trend data as well if it is not empty
    if ~isempty(trend_data)
        trend_data = trend_data(is_time);
    end

    % for returning
    cropped_dates = date_vector;

    % use datenums as x-axis
    x_values = datenum(date_vector);

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

    is_okay  = ~(isnan(data));

    % anything else but NaNs?
    if ( isempty(data(is_okay)) )
        anomalies = nan(size(data));
        return
    end

% ------------ %
% DETREND DATA %
% ------------ %

    if isempty(trend_data)

%          disp('calculating linear trend');

        % need to have NaN-free data for regression
        valid_data  = data(is_okay);
        valid_times = x_values(is_okay);

        % fitting: use only good values
        vals_poly = polyfit(valid_times,valid_data,1);

        % evaluation: use ALL values
        trend_data = polyval(vals_poly,x_values);

    % trend data compatible with data?
    else
        if ( length(data) ~= length(trend_data) )
            error('Data and trend data have different sizes.');
        end
    end

    % assign anomalies relative to the trend to
    % the correct piece of monthlyAnomalies
    detrended = data - trend_data;

% ------------------------------------ %
% REMOVE CLIMATOLOGICAL SEASONAL CYCLE %
% ------------------------------------ %

    anomalies = zeros(size(data));
    seasonal_cycle = zeros(1,12);

    for month = 1:12

        % timing
        is_month = ( date_vector(:,2) == month );


        % climatology
        climatology = nanmean(detrended(is_month));

        % assign anomalies
        anomalies(is_month)   = detrended(is_month)-climatology;
        seasonal_cycle(month) = nanmean(data(is_month));

    end