R/HW_index_3con.R

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

#' HWs of Three conditions
#'
#' @param ... other parameters to [HW_3con()]
#' @param date if NULL, `interval` will not use
#' 
#' @details
#' `T1 = quantile(x, p_high)`, `T2 = quantile(x, p_low)`.
#' 
#' HWs need to satisfy the following three conditions simultaneously:
#'   1. Everyday > T2;
#'   2. 3+ consecutive days > T1;
#'   3. Avg Tmax > T1 for whole period
#' 
#' @export
HW_3con <- function(x, p_high = 0.975, p_low = 0.81, mink =3,
    TRS_high = NULL, TRS_low = NULL,
    date = NULL,
    p1 = 1961,
    p2 = 1990,
    method = c("dt", "map"),
    diff = FALSE,
    ...)
{
    # if(is.matrix(x)) x <- heat_index(x[,1], x[,2])
    if (is.null(TRS_high)) TRS_high = cal_threshold(x, date, p_high)
    if (is.null(TRS_low)) TRS_low = cal_threshold(x, date, p_low)

    if (diff == TRUE) return(x - TRS_high) # might lead to error

    if (all(is.na(x))) return(rep(FALSE, length(x)))
    x[is.na(x)] <- -999L
    # val_max = max(x, na.rm = TRUE)
    # if (is.infinite(val_max) || val_max < TRS_high) return(rep(FALSE, length(x)))

    if (method[1] == "map") {
        x[x < TRS_low] <- NA
        grps <- split(x, cumsum(c(TRUE, diff(is.na(x)) != 0))) # 分组
        # HW_2con links to HW_2con.cpp
        status <- purrr::map_if(grps, ~ !all(is.na(.x)), ~ HW_2con(.x, TRS_high, mink)) %>%
            unlist() %>%
            as.logical()
        status[is.na(status)] <- FALSE
    } else if (method[1] == "dt") {
        if (is.null(date)) date = seq_along(x)
        dt <- data.table(date, x, status = FALSE)

        # HW_2con links to HW_2con.cpp
        dt[x > TRS_low,
            status := HW_2con(x, TRS_high, mink),
            by = cumsum(c(TRUE, diff(date) != 1L))]
        status = dt$status
    }
    status
}

#' @export
HW_3con_vec <- HW_3con

#' @export
#' @rdname HW_3con
HW_3con_rcpp <- function(x, p_high = 0.975, p_low = 0.81, mink = 3,
    TRS_high = NULL, TRS_low = NULL,
    date = NULL,
    p1 = 1961,
    p2 = 1990,
    # method = c("dt", "map"),
    # diff = FALSE,
    ...)
{
    # if (is.null(TRS_high)) TRS_high = quantile(x, p_high, na.rm = TRUE)
    # if (is.null(TRS_low)) TRS_low = quantile(x, p_low, na.rm = TRUE)
    # if(is.matrix(x)) x <- heat_index(x[,1], x[,2])
    if(all(is.na(x))) return(rep(FALSE, length(x)))
    if (is.null(TRS_high)) TRS_high = cal_threshold(x, date, p_high)
    if (is.null(TRS_low)) TRS_low = cal_threshold(x, date, p_low)
    x[is.na(x)] <- -999L
    # if (diff == TRUE) return(x - TRS_high) # might lead to error
    r = HW_episodes(x, TRS_high, TRS_low, mink = 3)
    rle_inverse(r$con3[1:3])
}

#' @param arr 3d array, with the dimension of `[nlon, nlat, ntime]`
#'
#' @rdname HW_3con
#' @export
HW_3con_array <- function(arr, ..., use.cpp = TRUE)
{
    FUN = ifelse(use.cpp, HW_3con_rcpp, HW_3con_vec)
    plyr::aaply(arr, c(1, 2), function(x) FUN(x, ...), .progress = "text")
}