R/HW_threshold.R

In CUG-hydro/heatwave: What the Package Does (One Line, Title Case)

# 分两步走
# 1. 先计算阈值;
# 2. 再计算算热浪指标

#' find dates in the interval of [p1, p2]
#'
#' @param p1 year_begin
#' @param p2 year_end
#' @param date Date object or character
#'
#' @return index vector
#' @export
find_interval <- function(date = NULL, p1 = 1961, p2 = 1990) {
    if (is.null(date)) return(seq_along(date))
    years = year(date)
    which(years >= p1 & years <= p2)
}

#' @export
cal_anomaly <- function(x, date = NULL, p1 = 1961, p2 = 1990) {
    ind = find_interval(date, p1, p2)
    x - mean(x[ind], na.rm = TRUE)
}

.probs <- c(75, 81, 90, 95, 97.5, 98, 99, 99.9, 99.99)/100
# mov: 考虑每天的阈值不同

#' HW threshold
#' 
#' @param x 3d array or numeric vector, with the dimension of `[lon, lat, ntime]`
#' @param method if `method = 'mov'``, 15 day moving average will be applied.
#' 
#' @importFrom rtrend movmean2_mat movmean
#' @export
cal_threshold <- function(x, date = NULL, probs = NULL, 
    p1 = 1961, p2 = 1990,
    method = c("none", "mov"),
    I_grid = NULL,
    halfwin = 7, ...)
{
    inds_ref = find_interval(date, p1, p2)
    ndim = length(dim(x))
    if (ndim == 3){
        mat = array_3dTo2d(x[, , inds_ref], I_grid)

        if (method[1] == "mov") {
            mat = movmean2_mat(mat, win_left = halfwin, win_right = halfwin)
        }
        if (is.null(probs)) probs = .probs
        TRS = rowQuantiles(mat, probs = probs, na.rm = TRUE)
        array_2dTo3d(TRS, I_grid, dim(x)[1:2])
    } else if (ndim == 0) {
        quantile(x[inds_ref], probs, na.rm = TRUE)
    }
}

#' @rdname cal_threshold
#' @export
cal_threshold_mov <- function(x, date = NULL, probs = NULL, 
    p1 = 1961, p2 = 1990,
    # method = c("fast", "low"),
    I_grid = NULL,
    halfwin = 7, ...)
{
    cal_threshold_doy <- function(doy) {
        ind_col = which(info$doy %in% find_adjacent_doy(doy, halfwin))
        rowQuantiles(mat[, ind_col], probs = probs, na.rm = TRUE) # TRS
    }

    if (is.null(probs)) probs = .probs
    inds_ref = find_interval(date, p1, p2)
    date = date[inds_ref]
    mat = if (length(dim(x)) == 3) array_3dTo2d(x[, , inds_ref], I_grid) else x[, inds_ref]

    info = date %>% data.table(date = ., date_md = format(., "%m-%d")) %>%
        mutate(doy = as.factor(date_md) %>% as.numeric())

    res = plyr::llply(1:366, cal_threshold_doy, .progress = "text")    
    map(1:length(probs), function(i){ map(res, ~.x[, i]) %>% do.call(cbind, .) }) %>% 
        set_names(paste0(probs*100, "%"))
}


find_adjacent_doy <- function(doy, halfwin = 7) {
    ind = (-halfwin:halfwin) + doy
    ind[ind > 366] %<>% subtract(366)
    ind[ind < 0]   %<>% add(366)
    ind
}


cal_threshold_low <- function(mat, dt, probs = NULL) {
    llply(1:nrow(mat), function(i) {
        dt$x = mat[i, ]
        thr <- slide(dt$x, ~.x, .before = 7, .after = 7) %>%
            set_names(dt$date_md) %>%
            utils::stack() %>%
            as.data.table() %>%
            .[, .(THR = quantile(values, probs, na.rm = TRUE)), by = ind]
        thr
    }, .progress = "text")
}