pixelTrend3D.m

function trend = pixelTrend3D(data)
%--------------------------------------------------------------------------
% pixelTrend3D computes the temporal trend for each pixel in a 3D array.
%--------------------------------------------------------------------------
% Input:
%   data: 3D numeric array (latitude x longitude x time)
% Output:
%   trend: 2D array (latitude x longitude) containing linear trends
%--------------------------------------------------------------------------
    
    % Validate input data
    if ndims(data) ~= 3
        error("Input data must be a 3D matrix (latitude x longitude x time)");
    end
    
    % Get the size of the data
    [rows, cols, ~] = size(data);

    % Initialize the trend array
    trend = NaN(rows, cols); % Use NaN to avoid uninitialized values

    % Calculate the trend for each (i, j) position
    for i = 1:rows
        for j = 1:cols
            % Extract the data for the current (i, j) position
            current_data = squeeze(data(i, j, :));
            
            % Check if sufficient non-NaN values exist for fitting
            valid_idx = ~isnan(current_data);
            
            if sum(valid_idx) > 1 % At least two valid points needed for polyfit
                coefficients = polyfit(find(valid_idx), current_data(valid_idx), 1);
                trend(i, j) = coefficients(1);
            end
        end
    end
end