R/Tinc_reach.R

In kongdd/CMIP5tools: Functions to Tidy and Analysis CMIP5/6 Data

# could improve a lot if wright in c++
# MAY HAVE BUG IF HAS NO HISTORICAL

#' Tinc_reach
#' 
#' @param Tinc `[ngrid, ntime]` increased Temperature
#' @param halfwin half size of moving window (year)
#' @param targets target warming levels
#' 
#' @return T_reach numeric vector `[ngrid]`
#' 
#' @keywords internal
#' @export
Tinc_reach <- function(Tinc, halfwin=10, targets = c(0.85, 1.5, 2, 3.5)) {
    tavg <- colMeans(Tinc)
    T_reach <- foreach(ij = 1:nrow(Tinc), .combine = rbind) %do% {
        y  <- Tinc[ij, ]
        xs <- phenofit::movmean(y, halfwin)
        
        tc <- foreach(Ti=targets, .combine = "c") %do% { which(xs >= Ti)[1] }
    } %>% set_colnames(targets)
    T_reach
}

#' @param Ti threshold of temperature 
#' 
#' @rdname Tinc_reach
#' @export
Tinc_reach.default <- function(Tinc, Ti){
    which(Tinc >= Ti)[1]
}

#' mean HW characteristics when reaching the target warming levels
#' 
#' @param mat matrix
#' @param T_reach Warming target reaching time.
#' @inheritParams Tinc_reach
#' 
#' @keywords internal
#' @export
mean_reach <- function(mat, T_reach, halfwin=10){
    nyear <- ncol(mat)
    # T_reach maybe NA value
    size <- dim(mat)
    # pos <- array(1:prod(size), dim = size)
    res <- foreach(Ti = iter(T_reach, "col"),.combine = cbind) %do% {
        temp <- foreach(t = Ti[,1], i = icount(),.combine = c, icount()) %do% {
            if (!is.na(t)) {
                i_start <- pmax(t - halfwin, 1)
                i_end   <- pmin(t + halfwin, nyear)
                # cbind(i, i_start, i_end)   
                vals = mat[i, i_start:i_end]
                # fix for FAR
                vals[is.infinite(vals)] = NA_real_
                mean(vals, na.rm = TRUE)
            } else {
                NA_real_
            }
        }
    }
    res %>% set_colnames(colnames(T_reach))
}

summary2 <- function(x){
    # x and Tinc should be smoothed by 20y-mov window first
    # I <- which(x >= Ti) %>% dplyr::first()
    x <- x[is.finite(x)]
    qs <- quantile(x, probs = c(0.25, 0.5, 0.75)) %>% 
        set_names(c("lower", "median", "upper")) %>% as.list()

    list(mean = mean(x), 
         max = max(x), 
         min = min(x), 
         sd = sd(x)) %>% 
    c(qs)
}