R/HW_index_song.R

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

#' At least k consecutive days > or ≥ T1, default is the former
#' 15-day moving sampling by year
#' @export
HW_1con_movingSPL_vec <- function(t, date, prob,
                                  k = 1, diff = FALSE,
                                  interval = 1961:1990, equate = FALSE) {
    if (is.matrix(t)) t <- map2_dbl(t[, 1], t[, 2], ~ heat_index(.x, .y))
    if (all(is.na(t)) == TRUE) {
        return(rep(FALSE, length(t)))
    }

    dt <- data.table(date, date_md = format(date, "%m-%d"), t)
    index <- dt[year(date) %in% interval, ]
    thr <- slide(index$t, ~.x, .before = 7, .after = 7) %>%
        set_names(index$date_md) %>%
        stack() %>%
        as.data.table() %>%
        .[, .(THR = quantile(values, prob, na.rm = TRUE)), by = ind]

    dt <- dt[thr, on = .(date_md = ind)][, status := "if"(equate, t >= THR, t > THR)]
    if (diff == TRUE) {
        return(dt[, t - THR])
    }
    
    if (k == 1) {
        return(replace_na(dt$status, FALSE))
    }

    run <- rle(dt$status)
    run$values[run$lengths < k] <- FALSE
    inverse.rle(run)
}

#' @export
HW_1con_movingSPL <- function(t_arr, rh_arr = NULL, ...) {
    if (!is.null(rh_arr)) t_arr <- abind(t_arr, rh_arr, along = 4)
    apply(t_arr, c(1, 2), function(x) HW_1con_movingSPL_vec(x, ...)) %>%
        aperm(c(2, 3, 1))
}

#' At least k consecutive days > or ≥ T1, default is the former
#' if percentile = True, THR = prob percentile of the data
#' input a arr, return a arr, including True or False
#'
#' background values is different
#'
#' @export
HW_1con_allSPL_vec <- function(t, date, k,
                               THR = NULL, prob = NULL,
                               diff = FALSE, interval = 1961:1990, equate = FALSE) {
    if (is.matrix(t)) t <- map2_dbl(t[, 1], t[, 2], ~ heat_index(.x, .y))
    if (all(is.na(t))) {
        return(rep(FALSE, length(t)))
    }

    dt <- data.table(date, t, status = FALSE)
    if (is.null(THR)) THR <- dt[year(date) %in% interval, quantile(t, prob, na.rm = TRUE)]
    if (diff == TRUE) {
        return(t - THR)
    }
    try(dt["if"(equate, t >= THR, t > THR),
        ":="(status = .N >= k),
        by = cumsum(c(TRUE, diff(date) != 1L))
    ], silent = TRUE)
    dt$status
}

## second method, maybe faster
# HW_1con_allSPL_vec <- function(t, date, k, THR = NULL, prob = NULL, interval = 1961:1990, equate = FALSE) {
#   if(all(is.na(t))) return(rep(FALSE, length(t)))
#
#   dt <- data.table(date, t)
#   if (!is.null(prob)) THR <- dt[year(date) %in% interval, quantile(t, prob, na.rm = TRUE)]
#   run <- rle('if'(equate, t >= THR, t > THR))
#   run$values[run$lengths < k] <- FALSE        #when k = 1, NA will not be replaced
#   inverse.rle(run)
# }
#' @export
HW_1con_allSPL <- function(t_arr, rh_arr = NULL, ...) {
    if (!is.null(rh_arr)) t_arr <- abind(t_arr, rh_arr, along = 4)
    apply(t_arr, MARGIN = c(1, 2), FUN = function(x) HW_1con_allSPL_vec(x, ...)) %>%
        aperm(c(2, 3, 1))
}