R/process_CMIP5.R

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

#' Process CMIP5 files by `FUN`
#' 
#' This function only for historical and RCP scenarios. This function is only 
#' used for temperature currently.
#' 
#' Historical scenario is forced to within 1850-2012
#' 
#' @inheritParams CMIP5Files_filter
#' @param lst_files list of CMIP5 files
#' @param outdir_TRS directory of Threshold RDS files
#' @param FUN see following examples
#' @param range `[lat_min, lat_max, lon_min, lon_max]`
#' @param outdir output directory
#' @param postfix postfix of output rds files
#' @param probs probabilities of non-exceeding
#' @param overwrite whether overwrite previous output
#' @param reverse whether process different scenarios in reverse order
#' @param .nchunk divide spatial grids into `.nchunk` pieces
#' @param ... other parameters to `FUN`
#' 
#' @examples
#' \dontrun{
#' FUN_rowMeans2 <- function(obj, ...) {
#'     apply_col(obj$value, year)
#' }
#' process_CMIP5(FUN, lst_files)
#' }
#' @export
process_CMIP5 <- function(
    FUN, 
    lst_files, outdir_TRS, outdir, 
    postfix = "_HW_index", 
    range, period = NULL, 
    varname = 'tasmax',
    value_range = NULL, probs, 
    overwrite = FALSE, 
    bigmemory = FALSE, 
    reverse = FALSE, 
    .nchunk = 6,
    ...) 
{
    check_dir(outdir)
    files_TRS <- dir(outdir_TRS, "*.RDS$", full.names = TRUE)
       
    r <- foreach(files = lst_files, scenario = names(lst_files), 
                 i=icount(), 
                 .combine = ) %do%{
        # setting period 
        time_info = get_period(scenario, period)

        ## MAIN SCRIPTS
        d_files <- CMIP5Files_filter(files, time_info$duration, time_info$period)
        lst_dfile <- d_files %>% split(., .$model)
        
        if (reverse) lst_dfile %<>% rev()
        # match TRS and ncfiles
        minfo <- match2(names(lst_dfile), 
                        get_model(files_TRS, prefix = "hold_", postfix = "\\."))
        
        foreach(d_file = lst_dfile[minfo$I_x],
                file_TRS = files_TRS[minfo$I_y], 
                .packages = c("CMIP5tools", "foreach", "data.table"),
                j = icount()) %do% {
            model  <- d_file$model[1]
            prefix <- sprintf('[i=%02d;j=%02d] %s, %s ', i, j, scenario, model)
            
            cat(sprintf("=========== %s =============\n", prefix))
            outfile <- sprintf("%s/%s_%s%s.RDS", outdir, scenario, model, postfix)
            
            if (file.exists(outfile) && !overwrite) return()

            obj_TRS <- readr::read_rds(file_TRS) # read TRS at here
            obj <- nc_merge(d_file, varname = varname, range, delta = 3, 
                            bigmemory = bigmemory, value_range = value_range, adjust_lon = TRUE)
            
            year <- obj$grid$date %>% str_year()
            nij  <- obj$grid$dim %>% {prod(.[1:2])}
            
            if (ncol(obj$data) != length(year)) {
                stop(sprintf("[e] date error!, %s\n", prefix))
            }
            
            ## write your FUNCTION at here
            cat(sprintf("APPLYING FUNCTION ... \n"))

            tryCatch({
                ## ERROR FOUND HERE
                r <- FUN(obj, obj_TRS, probs, prefix, .nchunk = .nchunk, ...)
                n_miss <- length(r$HW$year) - ncol(r$HW$duration$`0.9`)
                if (n_miss > 0) {
                    browser()
                }

                cat(sprintf("SAVING result ... \n"))
                saveRDS(r, outfile)
            }, error = function(e){
                message(sprintf("[e] %s: %s\n", prefix, e$message))
            }, finally = {
                rm(obj)
                gc(); gc()
                gc_cluster(.nchunk)
            })
        }
    }
}