scripts/ChinaHW_cluster/data02_Observed_TRS_HI&noHI.R

library(foreach)
library(nctools)
library(matrixStats)
InitCluster(6, kill = FALSE)

files_HI = dir("OUTPUT", "Observed_HI_V2_.*.nc", full.names = TRUE)
files_noHI = dir("data-raw", "CN05.1_T", full.names = TRUE)
names = paste(rep(c("HI", "noHI"), each = 3), c("Tavg", "Tmax", "Tmin"), sep = "_")
files = c(files_HI, files_noHI) %>% set_names(names)

I_grid = grid_d025$id %>% which.notna()
l_TRS = foreach(file = files_HI) %dopar% {
    l <- read_nc(file)
    dates = l$dates %>% as.Date.POSIXct()
    cal_threshold_mov(l$arr, dates, I_grid = I_grid)
}
save(l_TRS, file = "OUTPUT/Observed_HI&noHI_V2_TRS_15d-mov.rda")

## no mov
# I_grid <- grid_d025$id %>% which.notna()
# InitCluster(6, kill = FALSE)
l_TRS <- foreach(file = files_HI) %dopar% {
    l <- read_nc(file, year_end = 2016)
    dates <- l$dates %>% as.Date.POSIXct()
    cal_threshold(l$arr, dates, I_grid = I_grid)
}
save(l_TRS, file = "OUTPUT/Observed_TRS_HI&noHI_V2.rda")

## 3. check result (post processing)
files = dir("OUTPUT", "Observed_TRS_*", full.names = TRUE)
if (0) {
    load(files[1])
    l <- l_TRS %>% set_names(names[1:3])
    load(files[2])
    l_TRS[1:3] = l
    saveRDS(l_TRS, "OUTPUT/Observed_TRS_15d-mov_HI&noHI_V2.RDS", compress = FALSE)
    # diff = l_TRS$HI_Tavg$`75%` - l$HI_Tavg$`75%`

    load(files[3])
    l <- l_TRS %>% set_names(names[1:3])
    load(files[4])
    l_TRS[1:3] = l
    saveRDS(l_TRS, "OUTPUT/Observed_TRS_HI&noHI_V2.RDS", compress = FALSE)

}
if (0) {
    # check the result
    l <- readRDS("OUTPUT/Observed_TRS_15d-mov_HI&noHI_V2.RDS")
    load(files[1])
    diff = l_TRS[[1]]$`75%` - l$HI_Tavg$`75%`

    l <- readRDS("OUTPUT/Observed_TRS_HI&noHI_V2.RDS")
    load(files[3])
    diff = l_TRS[[1]] - l$HI_Tavg
    summary(as.numeric(diff))
}
CUG-hydro/heatwave documentation built on Dec. 17, 2021, 1:53 p.m.