R/process_raster2.R

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

#' combines historical and RCP scenario, then process
#' 
#' calculate certain warming reaching time and spatial average HW-index
#' 
#' @param type_TRSs can be some of `c("his30", "pi200")`
#' @param scenarios selected scenarios
#' @param FUN function to process `raster` object
#' @param year_begin begining year for processing
#' @param year_end ending year for processing
#' 
#' @param ... other parameters to `FUN`
#' 
#' @note This function is only for current file structure.
#' @export
process_raster2 <- function(
    type_TRSs = "pi200", 
    scenarios = c("historical", "RCP26", "RCP85"), # ,"RCP45"
    FUN = FUN_HW_reach, 
    year_begin = NULL, 
    year_end = 2100,
    ...)
{
    k <- 0
    type_TRSs %<>% set_names(., .)
    
    # different type_TRS
    res <- foreach(type_TRS = type_TRSs) %do% {
        lst_TRSs <- readRDS(sprintf("../DATA/HW_index/China_Threshold_%s.rds", type_TRS))
        
        # different variables
        foreach(varname = varnames) %do% {
            fprintf('[I] running [%s-%s] ... \n', type_TRS, varname)
            k <- k + 1
            # if (k <= 1) next()
            
            prefix <- sprintf('%s-%s', type_TRS, varname)
            files_HW <- dir(sprintf("../DATA/HW_index/China/%s",  prefix), "*.RDS", full.names = TRUE)
            lst_TRS <- lst_TRSs[[varname]]
            # lst_HW <- readRDS(sprintf("../DATA/HW_index/China_HW_index_%s_%s.rds",  type_TRS, varname))
            
            info <- gsub(prefix, "",  basename(files_HW)) %>% str_extract("(?<=_).*(?=_)") %>% {
                scenario <- str_extract(., ".*(?=_)")
                model <- str_extract(., "(?<=_).*")
               
                d <- data.table(file = files_HW, scenario, model = model, I_x = seq_along(files_HW))
                d$I_trs  <- match2(model, names(lst_TRS))$I_y 
                d %>% split(., .$scenario) %>% .[scenarios]
            }
            info_hist <- info$historical
            lst_info.rcp <- info[intersect(c("RCP26", "RCP45", "RCP85"), names(info))]

            # different RCP scenarios, combine historical by RCP
            temp <- foreach(info_rcp = lst_info.rcp, i = icount()) %do% {
                name_rcp <- names(lst_info.rcp)[i]
                minfo <- match2(info_hist$model, info_rcp$model) 
                
                foreach(j = 1:nrow(minfo)) %dopar% {
                    runningId(j, prefix = name_rcp)
                    r_hist <- read_rds(info_hist$file[minfo$I_x[j]])
                    r_rcp  <- read_rds(info_rcp$file[minfo$I_y[j]])
                    r <- combine_raster2(r_hist, r_rcp, year_begin = year_begin, year_end = year_end)
                    
                    I_trs <- info_hist$I_trs[minfo$I_x[j]]
                    TRS <- lst_TRS[[I_trs]]
                    # p <- profvis::profvis(d <- FUN(r, TRS, ...))
                    d <- FUN(r, TRS, ...)
                    d
                } %>% set_names(minfo$x)
            }
        }
    }
    return(res)
}

# Figure2 spatial
#' FUN of `process_raster2`
#' @name FUN_raster2
#' 
#' @description
#' - `FUN_HW_reach` : Reaching time of certain warming level
#' - `FUN_HW_annual`: annual time-series of regional mean
#' 
#' @param x raster2 obj of HW-index
#' @param TRS raster2 obj of TRS
#' 
#' @examples
#' \dontrun{
#' r <- FUN_HW_reach(x, TRS)
#' }
#' @keywords internal
NULL

#' @rdname FUN_raster2 
#' @export
FUN_HW_reach <- function(x, TRS){
    year <- colnames(x$T_annual) %>% as.numeric()
    Tinc <- x$T_annual - TRS$Tmean_PI
    T_reach <- Tinc_reach(Tinc)

    indices <- c("PR", "FAR") %>% set_names(., .)
    HW <- foreach(indice = indices) %do% {
        foreach(mat = x$HW[[indice]]) %do% {
            if (indice == "FAR") mat[!is.finite(mat)] <- -1 # missing_FAR
            mean_reach(mat, T_reach)
        }
    }
    HW$year <- year

    grid <- x[intersect(names(x), c("grid", "grid.origin"))]
    c(grid, listk(HW, T_reach)) %>% structure(class = "raster2")
}

# Figure1 temporal and Figure4 ECOF: temporal anormaly
# Interpolate and calculate regional mean
#' 
#' @rdname FUN_raster2 
#' @export
FUN_HW_annual <- function(x, TRS, gridInfo, range, p2 = NULL, 
    mutliple_regions = FALSE, 
    ...) 
{
    year <- colnames(x$T_annual) %>% as.numeric()
    x$HW$year <- year
    Tinc <- x$T_annual - TRS$Tmean_PI
    
    # show = TRUE, sp.layout =sp_arc_CH,
    Tinc.regional <- region_mean(Tinc, x$grid, gridInfo, range, 
        mutliple_regions = mutliple_regions, years = year)
    I <- match("value", colnames(Tinc.regional))
    colnames(Tinc.regional)[I] <- "Tinc"
    
    HW.regional <- region_mean_raster2(x$HW, x$grid, gridInfo, range, 
        mutliple_regions = mutliple_regions, p2 = p2)

    d <- merge(Tinc.regional, HW.regional)
    d
}