calcIndexAveraged.m 3.84 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
function index = calcIndexAveraged(monthlyAnomalies,latGrid,lonGrid,gridArea,latBoundaries,lonBoundaries)

% Calculates a climate index as area average over an area with given boundaries.
%
% Usage: index = calcIndexAveraged(monthlyAnomalies,latGrid,lonGrid,gridArea,latBoundaries,lonBoundaries)
%        index:            Time series of area averages (the "climate index").
%                          Same length as time dimension of monthlyAnomalies.
%        monthlyAnomalies: Raw data. Monthly anomalies.
%                          Dimensions: lon x lat x time
%        latGrid:          Latitude specification for monthlyAnomalies.
%                          Dimensions: lon x lat
%        lonGrid:          Same as latGrid, but for longitude.
%        gridArea:         Indicates the area of the individual grid boxes.
%                          Dimensions: lon x lat
%        latBoundaries:    Vector of length 2 that specifies the latitude boundaries of the area.
%                          latBoundaries(1): Minimum latitude
%                          latBoundaries(2): Maximum latitude
%        lonBoundaries:    Same as latBoundaries but for longitude.

% -------- %
% CHECKING %
% -------- %

    % Grid Sizes
    if ( ( size(monthlyAnomalies,1) ~= size(latGrid,1) ) | ...
         ( size(monthlyAnomalies,1) ~= size(lonGrid,1) ) | ...
         ( size(monthlyAnomalies,2) ~= size(latGrid,2) ) | ...
         ( size(monthlyAnomalies,2) ~= size(lonGrid,2) ) )
        error('Dimension between either of latGrid, lonGrid, and monthlyAnomalies!');
    end

    % Longitude format: grid
    if ( max(max(lonGrid)) > 180 )
        % disp('Reformatting: Maximum grid longitude > 180');
        lonGrid(lonGrid > 180) = lonGrid(lonGrid > 180) - 360;
    end

    % Longitude format: boundaries
    if ( max(lonBoundaries) > 180 )
        % disp('Reformatting: Maximum boundary longitude > 180');
        lonBoundaries(lonBoundaries > 180) = lonBoundaries(lonBoundaries > 180) - 360;
    end

% ------------------------------------------------ %
% FIND GRID POINTS CORRESPONDING TO THE GIVEN AREA %
% ------------------------------------------------ %

    % work with logical arrays the sime size as the data
    isLat = ( ( latGrid >= latBoundaries(1) ) & ( latGrid <= latBoundaries(2) ) );

    % does the longitudinal specification cross the date line?
    if ( ( lonBoundaries(1) > lonBoundaries(2) ) & ...
         ( sign(lonBoundaries(1)) == 1 ) & ( sign(lonBoundaries(2)) == -1 ) )
        % warning('Specified longitude boundaries cross the date line.');

        isLon_west = ( ( lonGrid >= lonBoundaries(1) ) & ( lonGrid <= 180 ) );
        isLon_east = ( ( lonGrid >= -180 )             & ( lonGrid <= lonBoundaries(2) ) );
        isLon      = ( isLon_west | isLon_east );

    % no crossing: first entry is smaller than the second one
    else
        isLon = ( ( lonGrid >= lonBoundaries(1) ) & ( lonGrid <= lonBoundaries(2) ) );
    end

% --------------- %
% CALCULATE INDEX %
% --------------- %

    % AREA MASK
        % initialize the mask:
        mask = ones(size(monthlyAnomalies));

        % multiply the (3D) mask with the logical array for the area:
        mask = bsxfun(@times,mask,(isLat & isLon));

        % set zeros to NaN:
        mask(mask == 0) = NaN;

    % mask anomalies:
    maskedAnomalies = monthlyAnomalies.*mask;

    % THE INDEX
        % total area of the region: Not a scalar but already of the same length as the
        % resulting time series --> can calculate the index time series
        % in one calculation step
        totalArea = ones(size(monthlyAnomalies,3),1).*sum(gridArea(isLat & isLon));

        % multiply anomalies with grid box area:
        weighted = bsxfun(@times,maskedAnomalies,gridArea);

        % THE INDEX:
        % area average: sum of the weighted anomalies divided by the total area
        index = squeeze(nansum(nansum(weighted,1),2))./totalArea;