%% Function cal_indice_GPM is used to calculate precipitation indices from raw precipitation data, which was established and tested on Matlab R2024b
%% For example, precipitation_data is a hourly time series (8760*1) of IMERG grid cell (19.95E,-40.05W) in 2001
%% Use the function to calculate precipitation characteristics of various intensities for each year and each grid cell
load('Precipitation series_example.mat');
precipitaiton_indice=cal_indice_GPM(precipitation_series, lon, 2001, 0);
%% Then you will get a matrix 'precipitaiton_indice' with dimension of 25*5*3
%% The first dimension means annual statistics (e.g., annual precipitation amount) and diurnal statistics from 00-23 local solar time
%% The second dimension means different precipitaiton types including total, drizzle, light, moderate, and heavy precipitation.
%% The third dimension means different metrics of precipitation including intensity, frequency, and amount.
%% See more details on our explanations below
function [precip_indice] = cal_indice_GPM(precip_period, lon, year, s)
%% 'precip_period' is a hourly precipitation time series derived from IMERG grid estimates (mm/h) in year (specified by parameter 'year')
%% 'lon' is the longitude of this grid cell
%% 's' defines the time scale for calculation, 0 for annual, 1 for DJF, 2 for MAM, 3 for JJA, and 4 for SON
    precip_threshold = [0.1, Inf; 0.1, 0.5; 0.5, 2.5; 2.5, 10; 10, Inf];
    precip_indice = ones(25, 5, 3) * nan;
    if s == 0
        days = yeardays(year);
    else
        switch s
            case 1  
                days = 31 + 28 + 31;
                if leapyear(year)
                    days = days + 1; 
                end
            case 2  
                days = 31 + 30 + 31;
            case 3  
                days = 30 + 31 + 31;
            case 4  
                days = 30 + 31 + 30;
            otherwise
                error('s must be integer between 0~4');
        end
    end
    valid_hours = sum(~isnan(precip_period));
    if valid_hours >= 0.8 * length(precip_period)
        for i = 1:5
            i0 = find(precip_period >= precip_threshold(i,1) & precip_period < precip_threshold(i,2));
            precip_indice(1,i,:) = [mean(precip_period(i0)), length(i0), sum(precip_period(i0))];
            precip_period1 = reshape(precip_period, 24, days);
            aa = ones(24, 3) * nan;
            for j = 1:24
                i1 = find(precip_period1(j,:) >= precip_threshold(i,1) & precip_period1(j,:) < precip_threshold(i,2));
                aa(j,:) = [mean(precip_period1(j,i1)), length(i1), sum(precip_period1(j,i1))];
            end
            k_local = (0:23) + round(lon / 15);
            k_local(k_local < 0) = k_local(k_local < 0) + 24;
            k_local(k_local > 23) = k_local(k_local > 23) - 24;
            precip_indice(2 + k_local, i, :) = aa;
        end
    end
end
