R/cal_Threshold.R

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

#' cal_Threshold
#' 
#' The quantile threshold and mean Temperature of piControl
#' see details in '/test/tasmax_s1_cal_threshold.R'.
#' 
#' @inheritParams threshold
#' 
#' @export
cal_Threshold <- function(d, range, outdir, 
    varname = 1L,
    bigmemory = FALSE, overwrite = FALSE, 
    .nchunk = 6)
{
    check_dir(outdir)
    outfile <- sprintf("%s/tasmax_Threshold_%s.RDS", outdir, d$model[1])
    # print(outfile)
    
    if (!file.exists(outfile) || overwrite) {
        tryCatch({
            # 1. constrain T in the range of [-100, 100]
            T_K_range <-  c(-100, 100) + 273.15            
            obj <- nc_merge(d, range, delta = 3, varname = varname, 
                            bigmemory = bigmemory, value_range = T_K_range)
            # # res_bilinear <- interp3d_bilinear(obj$dim, obj$data, range = range, cellsize_x = 1)
            probs   <- c(0.9, 0.95, 0.975, 0.99, 0.995, 0.9975, 0.999, 0.9995, 0.99975, 0.9999)
            mat_big <- obj$data
            
            lst_trs <- threshold(mat_big, probs, .nchunk, bigmemory = FALSE) %>% 
                set_names(c("TRS", "Tmean_PI"))

            ## PI control mean temperature
            r <- c(obj[1:2], lst_trs)
            
            cat(sprintf("writing %s\n", basename(outfile)))
            saveRDS(r, outfile)
            
            rm(obj); gc()
            r
        }, error = function(e){
            message(sprintf("[e] %s: %s", basename(outfile), e$message))
        })
    }
}

#' threshold of exceeding
#' 
#' @param mat_big temperature matrix, in the dimension of `[ngtid, ntime]`
#' @param .nchunk divide tasks into `.nchunk` pieces
#' @param probs probability of exceeding
#' 
#' @return 
#' * `TRS`: threshold
#' * `Tmean_ref`: mean temperature of reference period
#' 
#' @examples
#' \dontrun{
#' lst_trs <- threshold(mat_big, probs, .nchunk, bigmemory = FALSE)
#' }
#' 
#' @export
threshold <- function(mat_big, probs = c(0.99, 0.995), .nchunk = 6, bigmemory = FALSE) {
    if (bigmemory) xdesc <- bigmemory::describe(mat_big)

    cat(sprintf('running rowQuantiles ...\n'))
    # TEMP: could be improved
    # TRS <- rowQuantiles(mat_big, probs = probs, na.rm = TRUE)
    inds <- chunk(1:nrow(mat_big), .nchunk)
    TRS <- foreach(ind = inds, .combine = rbind, .packages = c("bigmemory")) %dopar% {
        if (bigmemory) mat_big <- bigmemory::attach.big.matrix(xdesc)
        apply(mat_big[ind, ], 1, quantile, probs = probs, na.rm = TRUE) %>% t()
    }
    
    cat(sprintf('running rowMeans2 ...\n'))
    # multi-year average
    Tmean_ref <- foreach(ind = inds, .combine = c) %dopar% {
        if (bigmemory) mat_big <- bigmemory::attach.big.matrix(xdesc)
        matrixStats::rowMeans2(mat_big[ind, ], na.rm = TRUE)
    }
    list(TRS = TRS, Tmean_ref = Tmean_ref)
}