diff --git a/R/QThreshold.R b/R/QThreshold.R index c4be966e4591a0087e1505cda1f6e837a42ac4d5..7d28b2a5a210be81f6a306054aad9efcd42228b7 100644 --- a/R/QThreshold.R +++ b/R/QThreshold.R @@ -3,6 +3,7 @@ #'Function to compute the corresponding quantile for a given value as required on function SU for the "Percentile_StressT" #'@param data A multidimensional array with named dimensions. #'@param threshold How many types of thresholds? When Experiments with several members are provided, threshold is not used. When Observations are provided it can be a single scalar. Should a user be able to provide an array? +#'@param ncores #'@importFrom s2dv Reorder #'@import multiApply #'@examples @@ -15,14 +16,14 @@ #'@export QThreshold <- function(data, threshold = NULL, memb_dim = 'member', - sdate_dim = 'sdate') { + sdate_dim = 'sdate', ncores = NULL) { dims <- dim(data) if (memb_dim %in% names(dims) && dims[which(names(dims) == memb_dim)] > 1) { res <- Apply(list(data), target_dims = c(memb_dim, sdate_dim), - fun = .qthreshold) + fun = .qthreshold, ncores = ncores) } else { res <- Apply(list(data), target_dims = sdate_dim, - fun = .qthreshold, threshold = threshold) + fun = .qthreshold, threshold = threshold, ncores = ncores) } return(res$output1) } diff --git a/R/SU.R b/R/SU.R index fd3cd722cfa897684b8d3f18c4e70278fd2cc6d2..f26b8f037a104f23b54b237eb6238d9b6997eff4 100644 --- a/R/SU.R +++ b/R/SU.R @@ -1,14 +1,17 @@ #' Index SU on multidimensional arrays #' SU35, SU36 and SU40 (Number of Heat Stress Days above a certain threshold) #' SU is the total count of days when the daily maximum temperatures exceed a certain threshold. - -### -# list of indices (threshold and period considered) -# SU35: 35C, seven months -# SU36: 36C, June 21 - September 21 -# SU40: 40C, June 21 - September 21 -# Spr32: 32C, April 21 - June 21 +#'@param data a multidimensional array with named dimensions +#'@param threshold a single scalar, vector or multidimensional array with named dimensions +#'@param time_dim a character string indicating the name of the temporal dimension to compute SU index +#'@param ncores a integer indicating the number of cores to use in parallel computation +#'@details list of indices (threshold and period considered) +#' SU35: 35C, seven months +#' SU36: 36C, June 21 - September 21 +#' SU40: 40C, June 21 - September 21 +#' Spr32: 32C, April 21 - June 21 ### +#'@import multiApply #'@examples #'tasmaxexp <- array(rnorm(120), #' c(member = 3, sdate = 4, ftime = 5, lat = 2)) @@ -16,29 +19,47 @@ #'tasmaxobs <- array(rnorm(40), #' c(member = 1, sdate = 4, ftime = 5, lat = 2)) #'threshold_obs <- QThreshold(tasmaxobs, 26) -#'res <- SU(tasmaxexp, tasmaxobs, threshold_exp, threshold_obs) +#'res <- SU(tasmaxexp, threshold_exp) +#'res <- SU(tasmaxobs, threshold_obs) +#'res <- SU(threshold_exp, threshold_obs) #'dim(res$su_exp) #'dim(res$su_obs) -SU <- function(exp, obs, threshold_exp, threshold_obs, threshold, time_dim = 'ftime') { +#'@export +SU <- function(data, threshold, time_dim = 'ftime', ncores = NULL) { #Checks missed - exp <- Apply(list(exp, threshold_exp), target_dims = list(time_dim, time_dim), - fun = .su)$output1 - obs <- Apply(list(obs, threshold_obs), target_dims = list(time_dim, time_dim), - fun = .su)$output1 - return(list(su_exp = exp, su_obs = obs)) -################################################# - # revising - - threshold_obs_reorder <- InsertDim(drop(threshold_obs), 1, dim(threshold_exp)['member']) - exp <- Apply(list(threshold_exp, threshold_obs_reorder), target_dims = list(time_dim, time_dim), - fun = .su)$output1 - - exp_nop <- Apply(list(exp), target_dims = time_dim, fun = .su, threshold = threshold)$output1 - obs <- Apply(list(obs), target_dims = time_dim, fun = .su, threshold = threshold)$output1 - - return(list(su_exp = exp, su_exp_nop = exp_nop, su_obs = obs)) -################################################# + if (is.null(dim(data))) { + dim(data) <- c(length(data)) + names(dim(data)) <- time_dim + } + if (is.null(dim(threshold)) && length(threshold) > 1) { + dim(threshold) <- c(length(threshold)) + names(dim(threshold)) <- time_dim + } + if (!is.character(time_dim)) { + stop("Parameter 'time_dim' must be a character string indicating the", + " name of the dimension to compute the index.") + } + if (!is.null(ncores)) { + if (!is.numeric(ncores)) { + ncores <- 1 + warning("Parameter 'ncores' must be numeric. Execution using 1 core.") + } + if (length(ncores) > 1) { + ncores <- ncores[1] + warning("Parameter 'ncores' has length greater than one and only", + " the first element will be used.") + } + } + threshold <- drop(threshold) + if (!is.null(dim(threshold)) && !all(dim(threshold) == 1)) { + data <- Apply(list(data, threshold), target_dims = list(time_dim, time_dim), + fun = .su, ncores = ncores)$output1 + } else { + data <- Apply(list(data), target_dims = list(time_dim), + fun = .su, threshold = threshold, ncores = ncores)$output1 + } + return(data) } -.su <- function(tasmax, threshold) { - result <- sum(tasmax > threshold) +.su <- function(data, threshold) { + result <- sum(data > threshold) } diff --git a/tests/testthat/test-SU.R b/tests/testthat/test-SU.R new file mode 100644 index 0000000000000000000000000000000000000000..eebfc122e7e43d43ee069cafc355fd61a48df9de --- /dev/null +++ b/tests/testthat/test-SU.R @@ -0,0 +1,22 @@ +context("Generic tests") +test_that("Sanity checks", { + data <- 1:20 + threshold <- 10 + expect_equal(SU(data, threshold), 10) + data <- array(1:40, c(x = 2, ftime = 20)) + expect_equal(SU(data, threshold), array(c(15, 15), c(x = 2))) + dim(threshold) <- c(member = 1, ftime = 1) + expect_equal(SU(data, threshold), array(c(15, 15), c(x = 2))) + expect_equal(SU(data, threshold, time_dim = 'x'), + array(c(rep(0,5), rep(2,15)), c(ftime = 20))) + expect_warning(SU(data, threshold, time_dim = 'x', ncores = 'Z'), + "Parameter 'ncores' must be numeric. Execution using 1 core.") + expect_warning(SU(data, threshold, time_dim = 'x', ncores = c(1,2)), + paste("Parameter 'ncores' has length greater than one and only", + "the first element will be used.")) + + expect_error(SU(data, threshold, time_dim = 3), + paste("Parameter 'time_dim' must be a character string", + "indicating the name of the dimension to compute the index.")) + +}) diff --git a/vignettes/BioClimaticIndicators.md b/vignettes/BioClimaticIndicators.md index 7a225ff735347aae7e28029edd3c48a5dff3a39f..2bb0fb98b0fcea260de23a70f4385ba87bb131be 100644 --- a/vignettes/BioClimaticIndicators.md +++ b/vignettes/BioClimaticIndicators.md @@ -7,81 +7,50 @@ This vignette shows how to compute Bioclimatic indicators developed under MEDGOL The number of Heat Stress Days is an index, called SU, reporting the total coung of days when daily maximum temperatures exceed 36 or 40°C between June 21st and September 21st. +Lines to be removed when the package is published: ``` # switch to ClimProjDiags local repo #setwd("/esarchive/scratch/nperez/git/ClimProjDiags") # My case library(s2dverification) +library(multiApply) +source("./R/SU.R") +source("./R/QThreshold.R") +``` -source('./R/SU.R') - +Lines to be removed when there is a call to load the data from esarchive: +``` pname <- '/esarchive/scratch/cchou/MEDGOLD/demo/output/su35/regrid2ERA5/data/Method_percentile/' load(paste0(pname, 'tasmax_era5grid_conservative_y9316.RData')) +``` +Data exploration: +``` dim(s5t_era5_9316) dim(era5t_era5_9316) tasmaxmod <- s5t_era5_9316[1:10, , 1:10, 1:3, 1:3] tasmaxobs <- era5t_era5_9316[, 1:10, 1:3, 1:3] -tasmaxobs_PStressT <- NULL stressT <- 35 -UsePercentile <- TRUE -ncores <- 2 -na.rm = TRUE - -library(multiApply) -data <- SU(tasmaxmod = tasmaxmod, tasmaxobs = tasmaxobs, tasmaxobs_PStressT = tasmaxobs_PStressT, - stressT = stressT, UsePercentile = UsePercentile, ncores = ncores, na.rm = na.rm) ``` -Exploring output: - +Explanation missed. ``` -names(data) -dim(data$sumod_WithPercentile) -dim(data$sumod_NoPercentile) -dim(data$suobs) -dim(data$tasmaxobs_PStressT) -dim(data$tasmaxobs_check) -dim(data$tasmaxmod_p) -dim(data$tasmaxmod_check_WithPercentile) -dim(data$tasmaxmod_check_NoPercentile) +threshold_obs <- QThreshold(tasmaxobs, stressT) +threshold_exp <- QThreshold(tasmaxmod) +su_exp <- SU(tasmaxmod, threshold_exp) # stressT instead of threshold_exp? +su_obs <- SU(tasmaxobs, threshold_obs) +su_thresholds <- SU(threshold_exp, threshold_obs) ``` -Result exploring output: +Exploring output: ``` -> names(data) -[1] "sumod_WithPercentile" "tasmaxobs_PStressT" -[3] "tasmaxmod_p" "tasmaxmod_check_WithPercentile" -[5] "sumod_NoPercentile" "suobs" -[7] "tasmaxmod_check_NoPercentile" "tasmaxobs_check" -> dim(data$sumod_WithPercentile) -member sdate lat lon - 10 24 3 3 -> dim(data$sumod_NoPercentile) -member sdate lat lon - 10 24 3 3 -> dim(data$suobs) -sdate lat lon - 24 3 3 -> dim(data$tasmaxobs_PStressT) -sdate ftime lat lon - 24 10 3 3 -> dim(data$tasmaxobs_check) -sdate ftime lat lon - 24 10 3 3 -> dim(data$tasmaxmod_p) -member sdate ftime lat lon - 10 24 10 3 3 -> dim(data$tasmaxmod_check_WithPercentile) -member sdate ftime lat lon - 10 24 10 3 3 -> dim(data$tasmaxmod_check_NoPercentile) -member sdate ftime lat lon - 10 24 10 3 3 +dim(su_exp) ``` +Visualization: + ``` -PlotEquiMap(data$sumod_WithPercentile[1,1,,], lon = 1:3, lat = 1:3, +PlotEquiMap(su_exp[1,1,,], lon = 1:3, lat = 1:3, filled.continents = FALSE) PlotEquiMap(data$sumod_NoPercentile[1,1,,], lon = 1:3, lat = 1:3, filled.continents = FALSE)