scripts/ChinaHW_cluster/data03_Observed_HW-cluster.R

In CUG-hydro/heatwave: What the Package Does (One Line, Title Case)

source("scripts/main_pkgs.R")
prj_dir = "/mnt/i/Research/cmip6/ChinaHW_cluster/"
load_all()
InitCluster(6, kill = FALSE)
setwd("../heatwave/")
{
    l_trs <- readRDS("OUTPUT/Observed_TRS_HI&noHI_V2.RDS")
    varnames <- c("Tavg", "Tmax", "Tmin")
    files_HI <- dir("OUTPUT", "Observed_HI_V2_.*.nc", full.names = TRUE) %>% set_names(varnames)
    files_noHI <- dir("data-raw", "CN05.1_T", full.names = TRUE) %>% set_names(varnames)

    names <- paste(rep(c("HI", "noHI"), each = 3), varnames, sep = "_")
    files <- c(files_HI, files_noHI) %>% set_names(names)

    I_grid <- grid_d025$id %>% which.notna()

    probs <- c(90, 95, 99, 99.9, 99.99) / 100 # Kong2020, JGR-A
    probs_name <- paste0(probs * 100, "%") %>% set_names(., .)

    date <- seq(as.Date("1961-01-01"), as.Date("2016-12-31"), by = "day")
}

## 1. 相对阈值法
if (0) {
    foreach(file = files, varname = c(varnames, varnames), prefix = names(files), k = icount()) %dopar% {
        if (k == 1) return()
        print(prefix)
        arr_raw <- read_nc(file, year_end = 2016)

        temp = foreach(prob_name = probs_name, prob = probs, i = icount()) %do% {
            outfile <- glue("OUTPUT/clusterId/dat-clusterId_{prefix}_perc{prob*100}.RDS")
            if (file.exists(outfile)) return()
            # arr_diff = arr_raw %-% TRS
            # arr_lgl <- default_false(arr_raw %>=% TRS)
            TRS <- l_trs[[prefix]][, , prob_name]
            system.time(clusterId <- HW_cluster(arr_lgl, date, ncell_connect = 16, TRS = TRS))
            saveRDS(clusterId, outfile, compress = TRUE)
        }
        gc(); gc()
    }
}

## 2. 滑动阈值法
if (0) {
    dim_xy <- c(283, 163)
    l_trs <- readRDS("OUTPUT/Observed_TRS_15d-mov_HI&noHI_V2.RDS")
    foreach(file = files, varname = c(varnames, varnames), prefix = names(files), k = icount()) %dopar% {
        # if (k == 1) return()
        print(prefix)
        arr_raw <- read_nc(file, year_end = 2016)

        temp = foreach(prob_name = probs_name, prob = probs, i = icount()) %do% {
            outfile <- glue("OUTPUT/clusterId/clusterId_movTRS_{prefix}_perc{prob*100}.RDS")
            if (file.exists(outfile)) return()
            TRS <- l_trs[[prefix]][[prob_name]] %>%
                array_2dTo3d(I_grid, dim_xy) # [[prefix]][[prob_name]][, , ]
            # arr_diff = arr_raw %-% TRS
            # arr_lgl <- default_false(arr_raw %>=% TRS)
            system.time(clusterId <- HW_cluster(arr_raw, date, ncell_connect = 16, TRS = TRS))
            saveRDS(clusterId, outfile, compress = TRUE)
        }
        gc(); gc()
    }
}

## 3. 绝对阈值法
if (0) {
    # dim_xy <- c(283, 163)
    # l_trs <- readRDS("OUTPUT/Observed_TRS_15d-mov_HI&noHI_V2.RDS")
    foreach(file = files, varname = c(varnames, varnames), prefix = names(files), k = icount()) %dopar% {
        # if (k == 1) return()
        print(prefix)
        if (!grepl("Tmax", prefix)) return()

        arr_raw <- read_nc(file, year_end = 2016)
        # temp = foreach(prob_name = probs_name, prob = probs, i = icount()) %do% {
        outfile <- glue("OUTPUT/clusterId/clusterId_absTRS(HW7)_{prefix}_perc{prob*100}.RDS")
        if (file.exists(outfile)) return()
        # TRS <- l_trs[[prefix]][, , prob_name]
        system.time(clusterId <- HW_cluster(arr_raw, date, ncell_connect = 16, TRS = 35))
        saveRDS(clusterId, outfile, compress = TRUE)
        # }
        gc(); gc()
    }
}

## 4. 双阈值法
# | HW10  | Tmax_2con_975 | K = 3  T1 = 97.5th percentile  T2 = 81th percentile |
# | HW11  | Tmax_2con_90  | K = 3  T1 = 90th percentile    T2 = 75th percentile |
# | HW12  | Tmin_2con_975 | K = 3  T1 = 97.5th percentile  T2 = 81th percentile |
# | HW13  | Tmin_2con_90  | K = 3  T1 = 90th percentile    T2 = 75th percentile |
if (1) {
    # l_trs <- readRDS("OUTPUT/Observed_TRS_HI&noHI_V2.RDS")
    foreach(file = files, varname = c(varnames, varnames), prefix = names(files), k = icount(6)) %do% {
        # if (k <= 3) return()
        print(prefix)
        if (grepl("Tmax", prefix)) return()
        arr_raw <- read_nc(file, year_end = 2016) # 3d array

        foreach(i = 1:2) %dopar% {
            if (i == 1) {
                p_low  = 0.75; p_high = 0.90
            } else {
                p_low  = 0.81; p_high = 0.975
            }
            outfile <- glue("OUTPUT/clusterId/clusterId_2conTRS(HW10)_{prefix}_perc{p_high*100}&{p_low*100}.RDS")
            if (file.exists(outfile)) next()

            print(outfile)
            arr_lgl <- HW_3con_array(arr_raw, date = date, p_low = p_low, p_high = p_high)
            clusterId <- HW_cluster(arr_lgl, date, ncell_connect = 16)
            saveRDS(clusterId, outfile, compress = TRUE)
        }
        gc(); gc()
    }
}