pixelstat2D.m

function [days_above_threshold, first_time_threshold, data_maxima] = pixelstat2D(data,threshold)
%------------------------------------------------------------------------------
% pixelstat2D computes pixelwise maxima of a data, 
% first day of occurrence of a threshold in each pixel
% and for how many days the data has remained above the thresold for each pixel
%------------------------------------------------------------------------------
% Input:
%   data: 3D numeric array (latitude x longitude x time)
%
% Output:
%   days_above_threshold: 
%   2D array (latitude x longitude) containing how many days 
%   the value was above the threshold (e.g. temperature above
%   5 degrees)
%
%   first_time_threshold: 
%   2D array (latitude x longitude) containing the
%   first time when the value rose abve the threshold
%
%   data_maxima: 
%   2D array (latitude x longitude) the maximum value of the
%   data in each of the locations
%------------------------------------------------------------------------------

    arguments
        data {mustBeNumeric, mustBe3D, mustBeWithin366} % should be a 3d input data (latitude x longitude x time)
        threshold (1,1) {mustBeNumeric} % can be negative or positive
    end

    % classify the data based on the matrix

    for i = 1:size(data,1)
        for j = 1:size(data,2)
            % Extract the time series for every pixel using the daily mean data
            timeseries = squeeze(data(i, j, :));
    
            % Count the number of days where data >= threshold
            days_above_threshold(i, j) = sum(timeseries >= threshold, 'omitnan');        
            
            % Find the first time index where data >= threshold
            idx = find(timeseries >= threshold, 1);
            
            % Store the time index if found, otherwise leave it as NaN
            if ~isempty(idx)
                first_time_threshold(i, j) = idx; % Store the time index
            end
            % find the maxima in each year for each of the pixel
            data_maxima(i, j) = max(data(i, j, :), [], 'omitnan');
        end
    end
end


% Custom validation function to check if the input data is 3D
function mustBe3D(data)
    if ndims(data) ~= 3
        error("Input data must be a 3D matrix (latitude x longitude x days)")
    end
end

function mustBeWithin366(data)
    if size(data,3) > 366
        error(['This is not a daily data for a year, The daily limit is 366, ...' ...
            'Please keep the data within the format (latitude x longitude x days)'])
    end
end