diff --git a/NAMESPACE b/NAMESPACE index 0a16d4d2dfb5888adf0007be1556694a7f87b899..7bed2f618f8cea4a56254e396088581a629f37c7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -7,6 +7,8 @@ export(CST_AccumulationExceedingThreshold) export(CST_MergeRefToExp) export(CST_PeriodAccumulation) export(CST_PeriodMean) +export(CST_PeriodPET) +export(CST_PeriodSPEI) export(CST_QThreshold) export(CST_SelectPeriodOnData) export(CST_Threshold) @@ -17,6 +19,8 @@ export(CST_WindPowerDensity) export(MergeRefToExp) export(PeriodAccumulation) export(PeriodMean) +export(PeriodPET) +export(PeriodSPEI) export(QThreshold) export(SelectPeriodOnData) export(SelectPeriodOnDates) @@ -25,7 +29,15 @@ export(TotalSpellTimeExceedingThreshold) export(TotalTimeExceedingThreshold) export(WindCapacityFactor) export(WindPowerDensity) +import(CSTools) +import(ClimProjDiags) +import(SPEI) +import(TLMoments) +import(lmom) +import(lmomco) +import(lubridate) import(multiApply) +import(zoo) importFrom(ClimProjDiags,Subset) importFrom(stats,approxfun) importFrom(stats,ecdf) diff --git a/R/PeriodPET.R b/R/PeriodPET.R new file mode 100644 index 0000000000000000000000000000000000000000..f365d39b444bca4aca8108e7b4f33523c896ba4c --- /dev/null +++ b/R/PeriodPET.R @@ -0,0 +1,409 @@ +#'Compute the Potential Evapotranspiration +#' +#'Compute the Potential evapotranspiration (PET) that is the amount of +#'evaporation and transpiration that would occur if a sufficient water source +#'were available. This function calculate PET according to the Thornthwaite, +#'Hargreaves or Hargreaves-modified equations. +#' +#'This function is build to work be compatible with other tools in +#'that work with 's2dv_cube' object class. The input data must be this object +#'class. If you don't work with 's2dv_cube', see PeriodPET. For more information +#'on the SPEI calculation, see functions CST_PeriodStandardization and +#'CST_PeriodAccumulation. +#' +#'@param data A named list with the needed \code{s2dv_cube} objects containing +#' the seasonal forecast experiment in the data element for each variable. +#' Specific variables are needed for each method used in computing the +#' Potential Evapotranspiration. See parameter 'pet_method'. The accepted +#' variables for each method are: for 'hargreaves', 'tmin' and 'tmax'; for +#' 'hargreaves_modified' are 'tmin', 'tmax' and 'pr'; for method 'thornthwaite' +#' 'tmean' is required. The units for temperature variables +#' ('tmin', 'tmax' and 'tmean') need to be in Celcius degrees; the units for +#' precipitation ('pr') need to be in mm/month. +#'@param dates An array of temporal dimensions containing the Dates of +#' 'data'. It must be of class 'Date' or 'POSIXct'. +#'@param lat A numeric vector containing the latitude values of 'data'. +#'@param pet_method A character string indicating the method used to compute +#' the potential evapotranspiration. The accepted methods are: +#' 'hargreaves' and 'hargreaves_modified', that require the data to have +#' variables tmin and tmax; and 'thornthwaite', that requires variable +#' 'tmean'. +#'@param time_dim A character string indicating the name of the temporal +#' dimension. By default, it is set to 'syear'. +#'@param leadtime_dim A character string indicating the name of the temporal +#' dimension. By default, it is set to 'time'. +#'@param lat_dim A character string indicating the name of the latitudinal +#' dimension. By default it is set by 'latitude'. +#'@param na.rm A logical value indicating whether NA values should be removed +#' from data. It is FALSE by default. +#'@param ncores An integer value indicating the number of cores to use in +#' parallel computation. +#' +#'@examples +#'dims <- c(time = 3, syear = 3, ensemble = 1, latitude = 1) +#' +#'exp_tasmax <- array(rnorm(360, 27.73, 5.26), dim = dims) +#'exp_tasmin <- array(rnorm(360, 14.83, 3.86), dim = dims) +#'exp_prlr <- array(rnorm(360, 21.19, 25.64), dim = dims) +#'end_year <- 2012 +#'dates_exp <- as.POSIXct(c(paste0(2010:end_year, "-08-16"), +#' paste0(2010:end_year, "-09-15"), +#' paste0(2010:end_year, "-10-16")), "UTC") +#'dim(dates_exp) <- c(syear = 3, time = 3) +#' +#'lat <- c(40) +#' +#'exp1 <- list('tmax' = exp_tasmax, 'tmin' = exp_tasmin, 'pr' = exp_prlr) +#' +#'res <- PeriodPET(data = exp1, lat = lat, dates = dates_exp) +#' +#'@import SPEI +#'@import lubridate +#'@import multiApply +#'@export +CST_PeriodPET <- function(data, pet_method = 'hargreaves', + time_dim = 'syear', leadtime_dim = 'time', + lat_dim = 'latitude', na.rm = FALSE, + ncores = NULL) { + # Check 's2dv_cube' + if (is.null(data)) { + stop("Parameter 'data' cannot be NULL.") + } + if (!all(sapply(data, function(x) inherits(x, 's2dv_cube')))) { + stop("Parameter 'data' must be a list of 's2dv_cube' class.") + } + # latitude + if (!any(names(data[[1]]$coords) %in% .KnownLatNames())) { + stop("Spatial coordinate names of parameter 'data' do not match any ", + "of the names accepted by the package.") + } + # Dates + dates_exp <- data[[1]]$attrs$Dates + if (!'Dates' %in% names(data[[1]]$attrs)) { + stop("Element 'Dates' is not found in 'attrs' list of 'data'. ", + "See 's2dv_cube' object description in README file for more ", + "information.") + } + lat_dim <- names(data[[1]]$coords)[[which(names(data[[1]]$coords) %in% .KnownLatNames())]] + + res <- PeriodPET(data = lapply(data, function(x) x$data), + dates = data[[1]]$attrs$Dates, + lat = data[[1]]$coords[[lat_dim]], + pet_method = pet_method, time_dim = time_dim, + leadtime_dim = leadtime_dim, lat_dim = lat_dim, + na.rm = na.rm, ncores = ncores) + # Add metadata + source_files <- lapply(data, function(x) {x$attrs$source_files}) + coords <- data[[1]]$coords + Dates <- data[[1]]$attrs$Dates + metadata <- data[[1]]$attrs$Variable$metadata + metadata_names <- intersect(names(dim(res)), names(metadata)) + suppressWarnings( + res <- CSTools::s2dv_cube(data = res, coords = coords, + varName = paste0('PET'), + metadata = metadata[metadata_names], + Dates = Dates, + source_files = source_files, + when = Sys.time()) + ) + return(res) +} + +#'Compute the Potential Evapotranspiration +#' +#'Compute the Potential evapotranspiration (PET) that is the amount of +#'evaporation and transpiration that would occur if a sufficient water source +#'were available. This function calculate PET according to the Thornthwaite, +#'Hargreaves or Hargreaves-modified equations. +#' +#'For more information on the SPEI calculation, see functions +#'PeriodStandardization and PeriodAccumulation. +#' +#'@param data A named list with the needed \code{s2dv_cube} objects containing +#' the seasonal forecast experiment in the data element for each variable. +#' Specific variables are needed for each method used in computing the +#' Potential Evapotranspiration. See parameter 'pet_method'. The accepted +#' variables for each method are: for 'hargreaves', 'tmin' and 'tmax'; for +#' 'hargreaves_modified' are 'tmin', 'tmax' and 'pr'; for method 'thornthwaite' +#' 'tmean' is required. The units for temperature variables +#' ('tmin', 'tmax' and 'tmean') need to be in Celcius degrees; the units for +#' precipitation ('pr') need to be in mm/month. +#'@param dates An array of temporal dimensions containing the Dates of +#' 'data'. It must be of class 'Date' or 'POSIXct'. +#'@param lat A numeric vector containing the latitude values of 'data'. +#'@param pet_method A character string indicating the method used to compute +#' the potential evapotranspiration. The accepted methods are: +#' 'hargreaves' and 'hargreaves_modified', that require the data to have +#' variables tmin and tmax; and 'thornthwaite', that requires variable +#' 'tmean'. +#'@param time_dim A character string indicating the name of the temporal +#' dimension. By default, it is set to 'syear'. +#'@param leadtime_dim A character string indicating the name of the temporal +#' dimension. By default, it is set to 'time'. +#'@param lat_dim A character string indicating the name of the latitudinal +#' dimension. By default it is set by 'latitude'. +#'@param na.rm A logical value indicating whether NA values should be removed +#' from data. It is FALSE by default. +#'@param ncores An integer value indicating the number of cores to use in +#' parallel computation. +#' +#'@examples +#'dims <- c(time = 3, syear = 3, ensemble = 1, latitude = 1) +#' +#'exp_tasmax <- array(rnorm(360, 27.73, 5.26), dim = dims) +#'exp_tasmin <- array(rnorm(360, 14.83, 3.86), dim = dims) +#'exp_prlr <- array(rnorm(360, 21.19, 25.64), dim = dims) +#'end_year <- 2012 +#'dates_exp <- as.POSIXct(c(paste0(2010:end_year, "-08-16"), +#' paste0(2010:end_year, "-09-15"), +#' paste0(2010:end_year, "-10-16")), "UTC") +#'dim(dates_exp) <- c(syear = 3, time = 3) +#' +#'lat <- c(40) +#' +#'exp1 <- list('tmax' = exp_tasmax, 'tmin' = exp_tasmin, 'pr' = exp_prlr) +#' +#'res <- PeriodPET(data = exp1, lat = lat, dates = dates_exp) +#' +#'@import SPEI +#'@import lubridate +#'@import multiApply +#'@export +PeriodPET <- function(data, dates, lat, pet_method = 'hargreaves', + time_dim = 'syear', leadtime_dim = 'time', + lat_dim = 'latitude', na.rm = FALSE, + ncores = NULL) { + + # Initial checks + ## data + if (!inherits(data, 'list')) { + stop("Parameter 'data' needs to be a named list with the needed variables.") + } + if (is.null(names(data))) { + stop("Parameter 'data' needs to be a named list with the variable names.") + } + if (any(sapply(data, function(x) is.null(names(dim(x)))))) { + stop("Parameter 'data' needs to be a list of arrays with dimension names.") + } + dims <- lapply(data, function(x) dim(x)) + first_dims <- dims[[1]] + all_equal <- all(sapply(dims[-1], function(x) identical(first_dims, x))) + if (!all_equal) { + stop("Parameter 'data' variables need to have the same dimensions.") + } + # lat + if (!is.numeric(lat)) { + stop("Parameter 'lat' must be numeric.") + } + if (!lat_dim %in% names(dims[[1]])) { + stop("Parameter 'data' must have 'lat_dim' dimension.") + } + if (any(sapply(dims, FUN = function(x) x[lat_dim] != length(lat)))) { + stop("Parameter 'lat' needs to have the same length of latitudinal", + "dimension of all the variables arrays in 'data'.") + } + + # data (2) + if (all(c('tmin', 'tmax', 'pr') %in% names(data))) { + # hargreaves modified: 'tmin', 'tmax', 'pr' and 'lat' + if (!(pet_method %in% c('hargreaves_modified', 'hargreaves'))) { + warning("Parameter 'pet_method' needs to be 'hargreaves' or ", + "'hargreaves_modified'. It is set to 'hargreaves_modified'.") + pet_method <- 'hargreaves_modified' + } + } else if (all(c('tmin', 'tmax') %in% names(data))) { + if (!(pet_method %in% c('hargreaves'))) { + warning("Parameter 'pet_method' will be set as 'hargreaves'.") + pet_method <- 'hargreaves' + } + } else if (c('tmean') %in% names(data)) { + # thornthwaite: 'tmean' (mean), 'lat' + if (!(pet_method == 'thornthwaite')) { + warning("Parameter 'pet_method' it is set to be 'thornthwaite'.") + pet_method <- 'thornthwaite' + } + } else { + stop("Parameter 'data' needs to be a named list with accepted ", + "variable names. See documentation.") + } + ## time_dim + if (!is.character(time_dim) | length(time_dim) != 1) { + stop("Parameter 'time_dim' must be a character string.") + } + if (!all(sapply(data, function(x) time_dim %in% names(dim(x))))) { + stop("Parameter 'time_dim' is not found in 'data' dimension.") + } + ## leadtime_dim + if (!is.character(leadtime_dim) | length(leadtime_dim) != 1) { + stop("Parameter 'leadtime_dim' must be a character string.") + } + if (!all(sapply(data, function(x) leadtime_dim %in% names(dim(x))))) { + stop("Parameter 'leadtime_dim' is not found in 'data' dimension.") + } + ## lat_dim + if (!is.character(lat_dim) | length(lat_dim) != 1) { + stop("Parameter 'lat_dim' must be a character string.") + } + if (!all(sapply(data, function(x) lat_dim %in% names(dim(x))))) { + stop("Parameter 'lat_dim' is not found in 'data' dimension.") + } + # dates + if (is.null(dates)) { + stop("Parameter 'dates' is missing, dates must be provided.") + } + if (!(is.Date(dates)) & !(is.POSIXct(dates))) { + stop("Parameter 'dates' is not of the correct class, ", + "only 'Date' and 'POSIXct' classes are accepted.") + } + if (!time_dim %in% names(dim(dates)) | !leadtime_dim %in% names(dim(dates))) { + stop("Parameter 'dates' must have 'time_dim' and 'leadtime_dim' ", + "dimension.") + } + if (!all(dim(data[[1]])[c(time_dim, leadtime_dim)] == + dim(dates)[c(time_dim, leadtime_dim)])) { + stop("Parameter 'dates' needs to have the same length as 'time_dim' ", + "and 'leadtime_dim' as 'data'.") + } + ## na.rm + if (!is.logical(na.rm) | length(na.rm) > 1) { + stop("Parameter 'na.rm' must be one logical value.") + } + ## ncores + if (!is.null(ncores)) { + if (!is.numeric(ncores) | any(ncores %% 1 != 0) | any(ncores < 0) | + length(ncores) > 1) { + stop("Parameter 'ncores' must be a positive integer.") + } + } + + # complete dates + dates_monthly <- .datesmask(dates) + + lat_mask <- array(lat, dim = c(1, length(lat))) + names(dim(lat_mask)) <- c('dat', lat_dim) + + # extract mask of NA locations to return to NA the final result + mask_na <- array(1, dim = dim(data[[1]])) + if (pet_method == 'hargreaves') { + varnames <- c('tmax', 'tmin') + mask_na[which(is.na(data$tmax))] <- 0 + mask_na[which(is.na(data$tmin))] <- 0 + } else if (pet_method == 'hargreaves_modified') { + varnames <- c('tmax', 'tmin', 'pr') + mask_na[which(is.na(data$tmax))] <- 0 + mask_na[which(is.na(data$tmin))] <- 0 + mask_na[which(is.na(data$pr))] <- 0 + } else if (pet_method == 'thornthwaite') { + varnames <- c('tmean') + mask_na[which(is.na(data$tmean))] <- 0 + } + + # replace NA with 0 + for (dd in 1:length(data)) { + data[[dd]][which(is.na(data[[dd]]))] <- 0 + } + + # prepare data + target_dims_data <- lapply(data[varnames], function(x) rep(c(leadtime_dim, time_dim), 1)) + pet <- Apply(data = c(list(lat_mask = lat_mask), data[varnames]), + target_dims = c(list(lat_mask = 'dat'), target_dims_data), + fun = .pet, + dates_monthly = dates_monthly, pet_method = pet_method, + leadtime_dim = leadtime_dim, time_dim = time_dim, + output_dims = c(leadtime_dim, time_dim), + ncores = ncores)$output1 + # reorder dims in pet_estimated + pos <- match(names(dim(data[[1]])), names(dim(pet))) + pet <- aperm(pet, pos) + + # restore original NAs from mask_na + pet[which(mask_na == 0)] <- NA + + return(pet) +} + +.pet <- function(lat_mask, data2, data3 = NULL, data4 = NULL, + dates_monthly, pet_method = 'hargreaves', + leadtime_dim = 'time', time_dim = 'syear') { + + dims <- dim(data2) + + # create a vector from data but adding 0 to achive complete time series + # of the considered period + # (starting in January of the first year) so that the solar radiation + # estimation is computed in each case for the correct month + + if (!is.null(data2)) { + data_tmp <- as.vector(data2) + data2 <- array(0, dim = length(dates_monthly)) + count <- 1 + for (dd in 1:length(dates_monthly)) { + if (dates_monthly[dd] == 1) { + data2[dd] <- data_tmp[count] + count <- count + 1 + } + } + rm(data_tmp) + } + if (!is.null(data3)) { + data_tmp <- as.vector(data3) + data3 <- array(0, dim = length(dates_monthly)) + count <- 1 + for (dd in 1:length(dates_monthly)) { + if (dates_monthly[dd] == 1) { + data3[dd] <- data_tmp[count] + count <- count + 1 + } + } + rm(data_tmp) + } + if (!is.null(data4)) { + data_tmp <- as.vector(data4) + data4 <- array(0, dim = length(dates_monthly)) + count <- 1 + for (dd in 1:length(dates_monthly)) { + if (dates_monthly[dd] == 1) { + data4[dd] <- data_tmp[count] + count <- count + 1 + } + } + rm(data_tmp) + } + if (pet_method == 'hargreaves') { + pet <- hargreaves(Tmin = as.vector(data3), Tmax = as.vector(data2), + lat = lat_mask, na.rm = FALSE, verbose = FALSE) + # line to return the vector to the size of the actual original data + pet <- array(pet[which(dates_monthly == 1)], dim = dims) + } + + if (pet_method == 'hargreaves_modified') { + pet <- hargreaves(Tmin = as.vector(data3), Tmax = as.vector(data2), + lat = lat_mask, Pre = as.vector(data4), na.rm = FALSE, + verbose = FALSE) + pet <- array(pet[which(dates_monthly == 1)], dim = dims) + } + + if (pet_method == 'thornthwaite') { + pet <- thornthwaite(as.vector(data2), lat = lat_mask, na.rm = TRUE, + verbose = FALSE) + # line to return the vector to the size of the actual original data + pet <- array(pet[which(dates_monthly == 1)], dim = dims) + } + return(pet) +} + + +.datesmask <- function(dates) { + ini <- as.Date(paste(lubridate::year(min(dates)), 01, 01, sep = '-')) + end <- as.Date(paste(lubridate::year(max(dates)), 12, 31, sep = '-')) + daily <- as.Date(ini:end) + monthly <- daily[which(lubridate::day(daily) == 1)] + dates_mask <- array(0, dim = length(monthly)) + for (dd in 1:length(dates)) { + ii <- which(monthly == as.Date(paste(lubridate::year(dates[dd]), + lubridate::month(dates[dd]), + 01, sep = '-'))) + dates_mask[ii] <- 1 + } + return(dates_mask) +} \ No newline at end of file diff --git a/R/PeriodSPEI.R b/R/PeriodSPEI.R new file mode 100644 index 0000000000000000000000000000000000000000..c2a04ecc0300c6f20ecd74ee0aa81c406421967d --- /dev/null +++ b/R/PeriodSPEI.R @@ -0,0 +1,1130 @@ +#'Compute the Standardised Precipitation-Evapotranspiration Index +#' +#'Calculation of the Standardised Precipitation-Evapotranspiration Index (SPEI) +#'that is a multiscalar drought index based on climatic data. It can be used for +#'determining the onset, duration and magnitude of drought conditions with +#'respect to normal conditions in a variety of natural and managed systems such +#'as crops, ecosystems, rivers, water resources, etc. The idea behind the SPEI +#'is to compare the highest possible evapotranspiration with the current water +#'availability. The SPEI uses the monthly (or weekly) difference between +#'precipitation and potential evapotranspiration. This represents a simple +#'climatic water balance which is calculated at different time scales to obtain +#'the SPEI. This function is build to work be compatible with other tools in +#'that work with 's2dv_cube' object class. The input data must be this object +#'class. If you don't work with 's2dv_cube', see PeriodSPEI. +#' +#'Next, some specifications for the calculation of this indicator will be +#'discussed. On the one hand, the model to be used to calculate potential +#'evapotranspiration is specified with the pet_method parameter (hargreaves, +#'hargraves modified or thornwhite). On the other hand, to choose the time scale +#'in which you want to accumulate the SPEI (SPEI3, SPEI6...) is done using the +#'accum parameter, where you must indicate the number of time steps you want to +#'accumulate throughout leadtime_dim. Since the accumulation is done for the +#'elapsed time steps, there will be no complete accumulations until reaching the +#'time instant equal to the value of the parameter. For this reason, in the +#'result, we will find that for the dimension where the accumulation has been +#'carried out, the values of the array will be NA since they do not include +#'complete accumulations. Also, there is a parameter to specify if the +#'standardization that can be TRUE or FALSE. If it is TRUE, the data is fit to a +#'probability distribution to transform the original values to standardized +#'units that are comparable in space and time and at different SPEI time scales. +#'The na.rm parameter is a logical parameter used to decide whether to remove +#'the NA values from the data before doing the calculation. It must be taken +#'into account that if na.rm == FALSE and there is some NA value in the specific +#'coordinates which the SPEI is computed, standardization cannot be carried out +#'for those coordinates and therefore, the result will be filled with NA for the +#'specific coordinates. However, when na.rm == TRUE, if the amount of data for +#'those specific coordinates is smaller than 4, it will not be possible to carry +#'out because we will not have enough data and the result will be also filled +#'with NAs for that coordinates. +#' +#'@param exp A named list with the needed \code{s2dv_cube} objects containing +#' the seasonal forecast experiment in the data element for each variable. +#' Specific variables are needed for each method used in computing the +#' Potential Evapotranspiration. See parameter 'pet_method'. The accepted +#' variables are: 'tmin' and 'tmax', required for methods 'hargreaves' and +#' 'hargreaves_modified' and 'tmean', required for method 'thornthwaite'. +#' Variable 'pr' is always needed. The units for temperature variables +#' ('tmin', 'tmax' and 'tmean') need to be in Celcius degrees; the units for +#' precipitation ('pr') need to be in mm/month. +#'@param exp_cor A named list with the needed \code{s2dv_cube} objects for each +#' variable in which the quantile PeriodSPEI should be applied. If it is not +#' specified, the PeriodSPEI is calculated from object 'exp'. +#'@param time_dim A character string indicating the name of the temporal +#' dimension. By default, it is set to 'syear'. +#'@param leadtime_dim A character string indicating the name of the temporal +#' dimension. By default, it is set to 'time'. +#'@param memb_dim A character string indicating the name of the dimension in +#' which the ensemble members are stored. When set it to NULL, threshold is +#' computed for individual members. +#'@param lat_dim A character string indicating the name of the latitudinal +#' dimension. By default it is set by 'latitude'. +#'@param accum An integer value indicating the number of months for the +#' accumulation for each variable. When it is greater than 1, the result will +#' be filled with NA until the accum time_dim dimension number due to the +#' accumulation to previous months. +#'@param ref_period A list with two numeric values with the starting and end +#' points of the reference period used for computing the index. The default +#' value is NULL indicating that the first and end values in data will be +#' used as starting and end points. +#'@param params A multidimensional array with named dimensions for computing the +#' SPEI. This option overrides computation of fitting parameters. It needs +#' to have the following dimensions: same leadtime dimension of 'exp' +#' (specified in 'leadtime_dim'); time dimension of length 1 (specified in +#' 'time_dim'); and a dimension named 'coef' with the length of the +#' coefficients needed for the used distribution (for 'Gamma' coef dimension is +#' of lenght 2, for 'log-Logistic' or 'PearsonIII' is of length 3). It can't +#' have member dimension (specified in 'memb_dim'). +#'@param pet_exp A multidimensional array containing the Potential +#' EvapoTranspiration data of 'exp'. It must have the same dimensions of the +#' variable 'pr' of 'exp'. If it is NULL it is calculated using the provided +#' variables with specified 'pet_method'. It is NULL by default. +#'@param pet_expcor A multidimensional array containing the Potential +#' EvapoTranspiration data of 'exp_cor'. It must have the same dimensions of +#' the variable 'pr' of 'exp_cor'. If it is NULL it is calculated using the +#' provided variables with specified 'pet_method'. It is NULL by default. +#'@param standardization A logical value indicating wether the standardization +#' is computed. +#'@param cross_validation A logical value indicating if cross validation is +#' done (TRUE), or not (FALSE). It only will be used when 'exp_cor' and +#' is not provided. It is FALSE by default. +#'@param pet_method A character string indicating the method used to compute +#' the potential evapotranspiration. The accepted methods are: +#' 'hargreaves' and 'hargreaves_modified', that require the data to have +#' variables tmin and tmax; and 'thornthwaite', that requires variable +#' 'tmean'. +#'@param method A character string indicating the standardization method used. +#' If can be: 'parametric' or 'non-parametric'. It is set to 'parametric' by +#' default. +#'@param distribution A character string indicating the name of the distribution +#' function to be used for computing the SPEI. The accepted names are: +#' 'log-Logistic' and 'Gamma'. It is set to 'log-Logistic' by default. The +#' 'Gamma' method only works when only precipitation is provided and other +#' variables are 0 because it is positive defined (SPI indicator). +#'@param param_error A numeric value with the error accepted. +#'@param handle_infinity A logical value wether to return Infinite values (TRUE) +#' or not (FALSE). +#'@param na.rm A logical value indicating whether NA values should be removed +#' from data. It is FALSE by default. If it is FALSE and there are NA values, +#' (if standardization is TRUE) all values of other dimensions except time_dim +#' and leadtime_dim will be set to NA directly. On the other hand, if it is +#' TRUE, if the data from other dimensions except time_dim and leadtime_dim is +#' not reaching 4 values, it is not enough values to estimate the parameters +#' and the result will include NA. +#'@param ncores An integer value indicating the number of cores to use in +#' parallel computation. +#' +#'@return An 's2dv_cube' object containing the SPEI multidimensional array in +#'element \code{data}. If 'exp_cor' is provided, only results from 'exp_cor' +#'will be provided. The parameters of the standardization will only be returned +#'if 'return_params' is TRUE. The SPEI will only be computed if +#''standardization' is TRUE. If 'standardization' is FALSE, only the climatic +#'water balance (precipitation minus evapotranspiration) will be returned. The +#'resultant arrays will have the same dimensions as the initial input data. The +#'other elements in the 's2dv_cube' will be updated with the combined +#'information of the input data arrays. +#' +#'@examples +#'dims <- c(var = 1, sday = 1, sweek = 1, syear = 6, time = 3, +#' latitude = 2, longitude = 1, ensemble = 25) +#'dimscor <- c(var = 1, sday = 1, sweek = 1, syear = 1, time = 3, +#' latitude = 2, longitude = 1, ensemble = 15) +#' +#'dates_exp <- as.POSIXct(c(paste0(2010:2015, "-08-16"), +#' paste0(2010:2015, "-09-15"), +#' paste0(2010:2015, "-10-16")), 'UTC') +#'dim(dates_exp) <- c(sday = 1, sweek = 1, syear = 6, time = 3) +#'dates_expcor <- as.POSIXct(c(paste0(2020, "-08-16"), paste0(2020, "-09-15"), +#' paste0(2020, "-10-16")), 'UTC') +#'dim(dates_expcor) <- c(sday = 1, sweek = 1, syear = 1, time = 3) +#'lat <- c(40,40.1) +#' +#'exp_tmax <- array(rnorm(900, 27.73, 5.26), dim = dims) +#'exp_tmin <- array(rnorm(900, 14.83, 3.86), dim = dims) +#'exp_pr <- array(rnorm(900, 21.19, 25.64), dim = dims) +#' +#'expcor_tmax <- array(rnorm(540, 29.03, 5.67), dim = dimscor) +#'expcor_tmin <- array(rnorm(540, 15.70, 4.40), dim = dimscor) +#'expcor_pr <- array(rnorm(540, 15.62, 21.38), dim = dimscor) +#' +#'exp <- list('tmax' = exp_tmax, 'tmin' = exp_tmin, 'pr' = exp_pr) +#'exp_cor <- list('tmax' = expcor_tmax, 'tmin' = expcor_tmin, +#' 'pr' = expcor_pr) +#' +#'exp <- lapply(exp, CSTools::s2dv_cube, coords = list(latitude = lat), +#' Dates = dates_exp) +#'exp_cor <- lapply(exp_cor, CSTools::s2dv_cube, coords = list(latitude = lat), +#' Dates = dates_expcor) +#' +#'res <- CST_PeriodSPEI(exp = exp, exp_cor = exp_cor) +#' +#'@import multiApply +#'@import ClimProjDiags +#'@import SPEI +#'@import zoo +#'@import TLMoments +#'@import lmomco +#'@import lmom +#'@import lubridate +#'@import CSTools +#'@export +CST_PeriodSPEI <- function(exp, exp_cor = NULL, + time_dim = 'syear', leadtime_dim = 'time', + memb_dim = 'ensemble', lat_dim = 'latitude', + accum = 1, ref_period = NULL, params = NULL, + pet_exp = NULL, pet_expcor = NULL, + standardization = TRUE, cross_validation = FALSE, + pet_method = 'hargreaves', method = 'parametric', + distribution = 'log-Logistic', + param_error = -9999, handle_infinity = FALSE, + return_params = FALSE, na.rm = FALSE, + ncores = NULL) { + + # Check 's2dv_cube' + if (is.null(exp)) { + stop("Parameter 'exp' cannot be NULL.") + } + if (!all(sapply(exp, function(x) inherits(x, 's2dv_cube')))) { + stop("Parameter 'exp' must be a list of 's2dv_cube' class.") + } + if (!is.null(exp_cor)) { + if (!all(sapply(exp_cor, function(x) inherits(x, 's2dv_cube')))) { + stop("Parameter 'exp_cor' must be a list of 's2dv_cube' class.") + } + } + # latitude + if (!any(names(exp[[1]]$coords) %in% .KnownLatNames())) { + stop("Spatial coordinate names of parameter 'exp' do not match any ", + "of the names accepted by the package.") + } + # Dates + dates_exp <- exp[[1]]$attrs$Dates + if (!'Dates' %in% names(exp[[1]]$attrs)) { + stop("Element 'Dates' is not found in 'attrs' list of 'exp'. ", + "See 's2dv_cube' object description in README file for more ", + "information.") + } + + if (!is.null(exp_cor)) { + if (!'Dates' %in% names(exp_cor[[1]]$attrs)) { + stop("Element 'Dates' is not found in 'attrs' list of 'exp_cor'. ", + "See 's2dv_cube' object description in README file for more ", + "information.") + } + } + + lat_dim <- names(exp[[1]]$coords)[[which(names(exp[[1]]$coords) %in% .KnownLatNames())]] + + res <- PeriodSPEI(exp = lapply(exp, function(x) x$data), + dates_exp = exp[[1]]$attrs$Dates, + lat = exp[[1]]$coords[[lat_dim]], + exp_cor = if (is.null(exp_cor)) {NULL} else { + lapply(exp_cor, function(x) x$data)}, + dates_expcor = exp_cor[[1]]$attrs$Dates, + time_dim = time_dim, leadtime_dim = leadtime_dim, + memb_dim = memb_dim, lat_dim = lat_dim, + accum = accum, ref_period = ref_period, params = params, + pet_exp = pet_exp, pet_expcor = pet_expcor, + standardization = standardization, + cross_validation = cross_validation, + pet_method = pet_method, method = method, + distribution = distribution, + param_error = param_error, handle_infinity = handle_infinity, + return_params = return_params, na.rm = na.rm, + ncores = ncores) + + if (!is.null(exp_cor)) { + source_files_expcor <- lapply(exp_cor, function(x) {x$attrs$source_files}) + source_files <- lapply(exp, function(x) {x$attrs$source_files}) + source_files <- c(exp = source_files, exp_cor = source_files_expcor) + coords <- exp_cor[[1]]$coords + Dates <- exp_cor[[1]]$attrs$Dates + metadata <- exp_cor[[1]]$attrs$Variable$metadata + metadata_names <- names(metadata) + } else { + source_files <- lapply(exp, function(x) {x$attrs$source_files}) + coords <- exp[[1]]$coords + Dates <- exp[[1]]$attrs$Dates + metadata <- exp[[1]]$attrs$Variable$metadata + metadata_names <- names(metadata) + } + + if (standardization) { + varname <- 'SPEI' + } else { + varname <- 'Precipitation minus accumulated PET' + } + + if (return_params & standardization) { + metadata_names <- intersect(names(dim(res[[1]])), metadata_names) + suppressWarnings( + res[[1]] <- CSTools::s2dv_cube(data = res[[1]], coords = coords, + varName = varname, + metadata = metadata[metadata_names], + Dates = Dates, + source_files = source_files, + when = Sys.time()) + ) + return(list(spei = res[[1]], params = res[[2]])) + } else { + metadata_names <- intersect(names(dim(res)), metadata_names) + suppressWarnings( + res <- CSTools::s2dv_cube(data = res, coords = coords, + varName = varname, + metadata = metadata[metadata_names], + Dates = Dates, + source_files = source_files, + when = Sys.time()) + ) + return(res) + } +} + + +#'Compute the Standardised Precipitation-Evapotranspiration Index +#' +#'Calculation of the Standardised Precipitation-Evapotranspiration Index (SPEI) +#'that is a multiscalar drought index based on climatic data. It can be used for +#'determining the onset, duration and magnitude of drought conditions with +#'respect to normal conditions in a variety of natural and managed systems such +#'as crops, ecosystems, rivers, water resources, etc. The idea behind the SPEI +#'is to compare the highest possible evapotranspiration with the current water +#'availability. The SPEI uses the monthly (or weekly) difference between +#'precipitation and potential evapotranspiration. This represents a simple +#'climatic water balance which is calculated at different time scales to obtain +#'the SPEI. +#' +#'Next, some specifications for the calculation of this indicator will be +#'discussed. On the one hand, the model to be used to calculate potential +#'evapotranspiration is specified with the pet_method parameter (hargreaves, +#'hargraves modified or thornwhite). On the other hand, to choose the time scale +#'in which you want to accumulate the SPEI (SPEI3, SPEI6...) is done using the +#'accum parameter, where you must indicate the number of time steps you want to +#'accumulate throughout leadtime_dim. Since the accumulation is done for the +#'elapsed time steps, there will be no complete accumulations until reaching the +#'time instant equal to the value of the parameter. For this reason, in the +#'result, we will find that for the dimension where the accumulation has been +#'carried out, the values of the array will be NA since they do not include +#'complete accumulations. Also, there is a parameter to specify if the +#'standardization that can be TRUE or FALSE. If it is TRUE, the data is fit to a +#'probability distribution to transform the original values to standardized +#'units that are comparable in space and time and at different SPEI time scales. +#'The na.rm parameter is a logical parameter used to decide whether to remove +#'the NA values from the data before doing the calculation. It must be taken +#'into account that if na.rm == FALSE and there is some NA value in the specific +#'coordinates which the SPEI is computed, standardization cannot be carried out +#'for those coordinates and therefore, the result will be filled with NA for the +#'specific coordinates. However, when na.rm == TRUE, if the amount of data for +#'those specific coordinates is smaller than 4, it will not be possible to carry +#'out because we will not have enough data and the result will be also filled +#'with NAs for that coordinates. +#' +#'@param exp A named list with multidimensional array objects containing +#' the seasonal forecast experiment in the data element for each variable. +#' Specific variables are needed for each method used in computing the +#' Potential Evapotranspiration. See parameter 'pet_method'. The accepted +#' variables are: 'tmin' and 'tmax', required for methods 'hargreaves' and +#' 'hargreaves_modified' and 'tmean', required for method 'thornthwaite'. +#' Variable 'pr' is always needed. The units for temperature variables +#' ('tmin', 'tmax' and 'tmean') need to be in Celcius degrees; the units for +#' precipitation ('pr') need to be in mm/month. +#'@param dates_exp An array of temporal dimensions containing the Dates of +#' 'exp'. It must be of class 'Date' or 'POSIXct'. +#'@param lat A numeric vector containing the latitude values of 'exp'. +#'@param exp_cor A named list with multidimensional array objects for each +#' variable in which the quantile PeriodSPEI should be applied. If it is not +#' specified, the PeriodSPEI is calculated from object 'exp'. +#'@param dates_expcor An array of temporal dimensions containing the Dates of +#' 'exp_cor'. It must be of class 'Date' or 'POSIXct'. +#'@param time_dim A character string indicating the name of the temporal +#' dimension. By default, it is set to 'syear'. +#'@param leadtime_dim A character string indicating the name of the temporal +#' dimension. By default, it is set to 'time'. +#'@param memb_dim A character string indicating the name of the dimension in +#' which the ensemble members are stored. When set it to NULL, threshold is +#' computed for individual members. +#'@param lat_dim A character string indicating the name of the latitudinal +#' dimension. By default it is set by 'latitude'. +#'@param accum accum An integer value indicating the number of months for the +#' accumulation for each variable. When it is greater than 1, the result will +#' be filled with NA until the accum time_dim dimension number due to the +#' accumulation to previous months. +#'@param ref_period A list with two numeric values with the starting and end +#' points of the reference period used for computing the index. The default +#' value is NULL indicating that the first and end values in data will be +#' used as starting and end points. +#'@param params A multidimensional array with named dimensions for computing the +#' SPEI. This option overrides computation of fitting parameters. It needs +#' to be of same time dimension (specified in 'time_dim') of 'exp' and a +#' dimension named 'coef' with the length of the coefficients needed for the +#' used distribution (for 'Gamma' coef dimension is of lenght 2, for +#' 'log-Logistic' or 'PearsonIII' is of length 3). It also needs to have a +#' leadtime dimension (specified in 'leadtime_dim') of length 1. +#'@param pet_exp A multidimensional array containing the Potential +#' EvapoTranspiration data of 'exp'. It must have the same dimensions of the +#' variable 'pr' of 'exp'. If it is NULL it is calculated using the provided +#' variables with specified 'pet_method'. It is NULL by default. +#'@param pet_expcor A multidimensional array containing the Potential +#' EvapoTranspiration data of 'exp_cor'. It must have the same dimensions of +#' the variable 'pr' of 'exp_cor'. If it is NULL it is calculated using the +#' provided variables with specified 'pet_method'. It is NULL by default. +#'@param standardization A logical value indicating wether the standardization +#' is computed. +#'@param cross_validation A logical value indicating if cross validation is +#' done (TRUE), or not (FALSE). It only will be used when 'exp_cor' and +#' is not provided. It is FALSE by default. +#'@param pet_method A character string indicating the method used to compute +#' the potential evapotranspiration. The accepted methods are: +#' 'hargreaves' and 'hargreaves_modified', that require the data to have +#' variables tmin and tmax; and 'thornthwaite', that requires variable +#' 'tmean'. +#'@param method A character string indicating the standardization method used. +#' If can be: 'parametric' or 'non-parametric'. +#'@param distribution A character string indicating the name of the distribution +#' function to be used for computing the SPEI. The accepted names are: +#' 'log-Logistic' and 'Gamma'. It is set to 'log-Logistic' by default. The +#' 'Gamma' method only works when only precipitation is provided and other +#' variables are 0 because it is positive defined (SPI indicator). +#'@param param_error A numeric value with the error accepted. +#'@param handle_infinity A logical value wether to return Infinite values (TRUE) +#' or not (FALSE). +#'@param return_params A logical value indicating wether to return parameters +#' array (TRUE) or not (FALSE). It is FALSE by default. +#'@param na.rm A logical value indicating whether NA values should be removed +#' from data. It is FALSE by default. If it is FALSE and there are NA values, +#' (if standardization is TRUE) all values of other dimensions except time_dim +#' and leadtime_dim will be set to NA directly. On the other hand, if it is +#' TRUE, if the data from other dimensions except time_dim and leadtime_dim is +#' not reaching 4 values, it is not enough values to estimate the parameters +#' and the result will include NA. +#'@param ncores An integer value indicating the number of cores to use in +#' parallel computation. +#' +#'@return An 's2dv_cube' object containing the SPEI multidimensional array in +#'element \code{data}. If 'exp_cor' is provided, only results from 'exp_cor' +#'will be provided. The parameters of the standardization will only be returned +#'if 'return_params' is TRUE. The SPEI will only be computed if +#''standardization' is TRUE. If 'standardization' is FALSE, only the climatic +#'water balance (precipitation minus evapotranspiration) will be returned. The +#'resultant arrays will have the same dimensions as the initial input data. +#' +#'@examples +#'dims <- c(var = 1, sday = 1, sweek = 1, syear = 6, time = 3, +#' latitude = 2, longitude = 1, ensemble = 25) +#'dimscor <- c(var = 1, sday = 1, sweek = 1, syear = 1, time = 3, +#' latitude = 2, longitude = 1, ensemble = 15) +#'exp_tmax <- array(rnorm(900, 27.73, 5.26), dim = dims) +#'exp_tmin <- array(rnorm(900, 14.83, 3.86), dim = dims) +#'exp_pr <- array(rnorm(900, 21.19, 25.64), dim = dims) +#' +#'expcor_tmax <- array(rnorm(540, 29.03, 5.67), dim = dimscor) +#'expcor_tmin <- array(rnorm(540, 15.70, 4.40), dim = dimscor) +#'expcor_pr <- array(rnorm(540, 15.62, 21.38), dim = dimscor) +#' +#'dates_exp <- as.Date(c(paste0(2010:2015, "-08-16"), +#' paste0(2010:2015, "-09-15"), +#' paste0(2010:2015, "-10-16"))) +#'dim(dates_exp) <- c(sday = 1, sweek = 1, syear = 6, time = 3) +#'dates_expcor <- as.POSIXct(c(paste0(2020, "-08-16"), +#' paste0(2020, "-09-15"), +#' paste0(2020, "-10-16")), 'UTC') +#'dim(dates_expcor) <- c(sday = 1, sweek = 1, syear = 1, time = 3) +#'lat <- c(40,40.1) +#' +#'exp <- list('tmax' = exp_tmax, 'tmin' = exp_tmin, 'pr' = exp_pr) +#'exp_cor <- list('tmax' = expcor_tmax, 'tmin' = expcor_tmin, +#' 'pr' = expcor_pr) +#'res <- PeriodSPEI(exp = exp, exp_cor = exp_cor, lat = lat, +#' dates_exp = dates_exp, dates_expcor = dates_expcor) +#' +#'@import multiApply +#'@import ClimProjDiags +#'@import SPEI +#'@import zoo +#'@import TLMoments +#'@import lmomco +#'@import lmom +#'@import lubridate +#'@export +PeriodSPEI <- function(exp, dates_exp, lat, + exp_cor = NULL, dates_expcor = NULL, + time_dim = 'syear', leadtime_dim = 'time', + memb_dim = 'ensemble', lat_dim = 'latitude', + accum = 1, ref_period = NULL, params = NULL, + pet_exp = NULL, pet_expcor = NULL, + standardization = TRUE, cross_validation = FALSE, + pet_method = 'hargreaves', method = 'parametric', + distribution = 'log-Logistic', + param_error = -9999, handle_infinity = FALSE, + return_params = FALSE, na.rm = FALSE, ncores = NULL) { + + # Initial checks + ## exp + if (!inherits(exp, 'list')) { + stop("Parameter 'exp' needs to be a named list with the needed variables.") + } + if (is.null(names(exp))) { + stop("Parameter 'exp' needs to be a named list with the variable names.") + } + if (any(sapply(exp, function(x) is.null(names(dim(x)))))) { + stop("Parameter 'exp' needs to be a list of arrays with dimension names.") + } + dims <- lapply(exp, function(x) dim(x)) + first_dims <- dims[[1]] + all_equal <- all(sapply(dims[-1], function(x) identical(first_dims, x))) + if (!all_equal) { + stop("Parameter 'exp' variables need to have the same dimensions.") + } + + ## exp_cor + if (!is.null(exp_cor)) { + if (!inherits(exp_cor, 'list')) { + stop("Parameter 'exp_cor' needs to be a named list with the needed ", + "variables if it is not NULL.") + } + if (is.null(names(exp_cor))) { + stop("Parameter 'exp_cor' needs to be a named list with the variable names.") + } + if (any(sapply(exp_cor, function(x) is.null(names(dim(x)))))) { + stop("Parameter 'exp_cor' needs to be a list of arrays with dimension names.") + } + dimscor <- lapply(exp_cor, function(x) dim(x)) + first_dimscor <- dimscor[[1]] + all_equal <- all(sapply(dimscor[-1], function(x) identical(first_dimscor, x))) + if (!all_equal) { + stop("Parameter 'exp_cor' variables need to have the same dimensions.") + } + } + # Variable checks + ## exp (2) + pet <- vector("list", 2) + if (!('pr' %in% names(exp))) { + stop("Variable 'pr' is not included in 'exp'.") + } + ## exp_cor (2) + if (!is.null(exp_cor)) { + if (!('pr' %in% names(exp_cor))) { + stop("Variable 'pr' is not included in 'exp_cor'.") + } + if (length(pet_method) == 1) { + pet_method <- rep(pet_method, 2) + } + } + + ## pet_exp + if (!is.null(pet_exp)) { + if (length(dim(exp[['pr']])) != length(dim(pet_exp))) { + stop("Parameter 'pet_exp' must have the same length of all the ", + "dimensions as variable 'pr' in 'exp'.") + } + if (!all(dim(exp[['pr']]) %in% dim(pet_exp))) { + stop("Parameter 'pet_exp' must have the same length of all the ", + "dimensions as variable 'pr' in 'exp'.") + } + if (any(names(dim(exp[['pr']])) != names(dim(pet_exp)))) { + pos <- match(names(dim(exp[['pr']])), names(dim(pet_exp))) + pet_exp <- aperm(pet_exp, pos) + } + pet[[1]] <- pet_exp + } else if (is.null(dates_exp)) { + stop("Parameter 'dates_exp' must be provided.") + } + ## pet_expcor + if (!is.null(exp_cor)) { + if (!is.null(pet_expcor)) { + if (length(dim(exp_cor[['pr']])) != length(dim(pet_expcor))) { + stop("Parameter 'pet_expcor' must have the same length of all the ", + "dimensions as variable 'pr' in 'exp_cor'.") + } + if (!all(dim(exp_cor[['pr']]) %in% dim(pet_expcor))) { + stop("Parameter 'pet_expcor' must have the same length of all the ", + "dimensions as variable 'pr' in 'exp_cor'.") + } + if (any(names(dim(exp_cor[['pr']])) != names(dim(pet_expcor)))) { + pos <- match(names(dim(exp_cor[['pr']])), names(dim(pet_expcor))) + pet_expcor <- aperm(pet_expcor, pos) + } + pet[[2]] <- pet_expcor + } else if (is.null(dates_expcor)) { + stop("Parameter 'dates_expcor' must be provided.") + } + } + + ## time_dim + if (!is.character(time_dim) | length(time_dim) != 1) { + stop("Parameter 'time_dim' must be a character string.") + } + if (!all(sapply(exp, function(x) time_dim %in% names(dim(x))))) { + stop("Parameter 'time_dim' is not found in 'exp' dimension.") + } + if (!is.null(exp_cor)) { + if (!all(sapply(exp_cor, function(x) time_dim %in% names(dim(x))))) { + stop("Parameter 'time_dim' is not found in 'exp_cor' dimension.") + } + } + ## leadtime_dim + if (!is.character(leadtime_dim) | length(leadtime_dim) != 1) { + stop("Parameter 'leadtime_dim' must be a character string.") + } + if (!all(sapply(exp, function(x) leadtime_dim %in% names(dim(x))))) { + stop("Parameter 'leadtime_dim' is not found in 'exp' dimension.") + } + if (!is.null(exp_cor)) { + if (!all(sapply(exp_cor, function(x) leadtime_dim %in% names(dim(x))))) { + stop("Parameter 'leadtime_dim' is not found in 'exp_cor' dimension.") + } + } + ## memb_dim + if (!is.character(memb_dim) | length(memb_dim) != 1) { + stop("Parameter 'memb_dim' must be a character string.") + } + if (!all(sapply(exp, function(x) memb_dim %in% names(dim(x))))) { + stop("Parameter 'memb_dim' is not found in 'exp' dimension.") + } + if (!is.null(exp_cor)) { + if (!all(sapply(exp_cor, function(x) memb_dim %in% names(dim(x))))) { + stop("Parameter 'memb_dim' is not found in 'exp_cor' dimension.") + } + } + + ## lat_dim + if (!is.character(lat_dim) | length(lat_dim) != 1) { + stop("Parameter 'lat_dim' must be a character string.") + } + if (!all(sapply(exp, function(x) lat_dim %in% names(dim(x))))) { + stop("Parameter 'lat_dim' is not found in 'exp' dimension.") + } + if (!is.null(exp_cor)) { + if (!all(sapply(exp_cor, function(x) lat_dim %in% names(dim(x))))) { + stop("Parameter 'lat_dim' is not found in 'exp_cor' dimension.") + } + } + + ## dates + if (is.null(dates_exp)) { + stop("Parameter 'dates_exp' is missing, dates must be provided.") + } + if (!(is.Date(dates_exp)) & !(is.POSIXct(dates_exp))) { + stop("Parameter 'dates_exp' is not of the correct class, ", + "only 'Date' and 'POSIXct' classes are accepted.") + } + if (!time_dim %in% names(dim(dates_exp)) | !leadtime_dim %in% names(dim(dates_exp))) { + stop("Parameter 'dates_exp' must have 'time_dim' and 'leadtime_dim' ", + "dimension.") + } + if (!all(dim(exp[[1]])[c(time_dim, leadtime_dim)] == + dim(dates_exp)[c(time_dim, leadtime_dim)])) { + stop("Parameter 'dates_exp' needs to have the same length as 'time_dim' ", + "and 'leadtime_dim' as 'exp'.") + } + + if (!is.null(exp_cor)) { + if (is.null(dates_expcor)) { + stop("Parameter 'dates_expcor' is missing, dates for 'exp_cor' must be ", + "provided if exp_cor is not NULL.") + } + if (!(is.Date(dates_expcor)) & !(is.POSIXct(dates_expcor))) { + stop("Element 'Dates' in 'exp_cor' is not of the correct class, ", + "only 'Date' and 'POSIXct' classes are accepted.") + } + if (!time_dim %in% names(dim(dates_expcor)) | !leadtime_dim %in% names(dim(dates_expcor))) { + stop("Parameter 'dates_expcor' must have 'time_dim' and 'leadtime_dim' ", + "dimension.") + } + if (!all(dim(exp_cor[[1]])[c(time_dim, leadtime_dim)] == + dim(dates_expcor)[c(time_dim, leadtime_dim)])) { + stop("Parameter 'dates_expcor' needs to have the same length as ", + "'time_dim' and 'leadtime_dim' as 'exp_cor'.") + } + } + + ## accum + if (accum > dim(exp[[1]])[leadtime_dim]) { + stop(paste0("Cannot compute accumulation of ", accum, " months because ", + "loaded data has only ", dim(exp[[1]])[leadtime_dim], " months.")) + } + + ## ref_period + if (!is.null(ref_period)) { + if (length(ref_period) != 2) { + warning("Parameter 'ref_period' must be of length two indicating the ", + "first and end years of the reference period. It will not ", + "be used.") + ref_period <- NULL + } else if (!all(sapply(ref_period, is.numeric))) { + warning("Parameter 'ref_period' must be a numeric vector indicating the ", + "'start' and 'end' years of the reference period. It will not ", + "be used.") + ref_period <- NULL + } else if (ref_period[[1]] > ref_period[[2]]) { + warning("In parameter 'ref_period' 'start' cannot be after 'end'. It ", + "will not be used.") + ref_period <- NULL + } else if (!all(unlist(ref_period) %in% year(dates_exp))) { + warning("Parameter 'ref_period' contain years outside the dates. ", + "It will not be used.") + ref_period <- NULL + } else { + years <- year(ClimProjDiags::Subset(dates_exp, along = leadtime_dim, + indices = 1)) + ref_period[[1]] <- which(ref_period[[1]] == years) + ref_period[[2]] <- which(ref_period[[2]] == years) + } + } + + ## standardization + if (!is.logical(standardization)) { + stop("Parameter 'standardization' must be a logical value.") + } + + ## param_error + if (!is.numeric(param_error)) { + stop("Parameter 'param_error' must be a numeric value.") + } + + ## handle_infinity + if (!is.logical(handle_infinity)) { + stop("Parameter 'handle_infinity' must be a logical value.") + } + + ## cross_validation + if (!is.logical(cross_validation)) { + stop("Parameter 'cross_validation' must be a logical value.") + } + if (!is.null(exp_cor)) { + if (cross_validation & standardization) { + warning("Detected 'cross_validation' = TRUE. It will be set as FALSE ", + "since 'exp_cor' is provided.") + cross_validation <- FALSE + } + } + + ## method + if (!(method %in% c('parametric', 'non-parametric'))) { + stop("Parameter 'method' must be a character string containing one of ", + "the following methods: 'parametric' or 'non-parametric'.") + } + + ## distribution + if (!(distribution %in% c('log-Logistic', 'Gamma', 'PearsonIII'))) { + stop("Parameter 'distribution' must be a character string containing one ", + "of the following distributions: 'log-Logistic', 'Gamma' or ", + "'PearsonIII'.") + } + ## ncores + if (!is.null(ncores)) { + if (!is.numeric(ncores) | any(ncores %% 1 != 0) | any(ncores < 0) | + length(ncores) > 1) { + stop("Parameter 'ncores' must be a positive integer.") + } + } + + ## na.rm + if (!is.logical(na.rm)) { + stop("Parameter 'na.rm' must be logical.") + } + + ## params + if (!is.null(params)) { + if (!is.numeric(params)) { + stop("Parameter 'params' must be numeric.") + } + if (!all(c(time_dim, leadtime_dim, 'coef') %in% names(dim(params)))) { + stop("Parameter 'params' must be a multidimensional array with named ", + "dimensions: '", time_dim, "', '", leadtime_dim, "' and 'coef'.") + } + if (distribution == "Gamma") { + if (dim(params)['coef'] != 2) { + stop("For '", distribution, "' distribution, params array should have ", + "'coef' dimension of length 2.") + } + } else { + if (dim(params)['coef'] != 3) { + stop("For '", distribution, "' distribution, params array should have ", + "'coef' dimension of length 3.") + } + } + } + + ## return_params + if (!is.logical(return_params)) { + stop("Parameter 'return_params' must be logical.") + } + + # Complete dates + dates <- .return2list(dates_exp, dates_expcor) + + # Compute PeriodSPEI + k = 0 + spei_res <- NULL + computed_pet <- FALSE + + for (data in .return2list(exp, exp_cor)) { + k = k + 1 + # Evapotranspiration estimation (unless pet is already provided) + if (is.null(pet[[k]]) | computed_pet) { + pet[[k]] <- PeriodPET(data = data, dates = dates[[k]], + lat = lat, pet_method = pet_method[k], + time_dim = time_dim, leadtime_dim = leadtime_dim, + lat_dim = lat_dim, na.rm = na.rm, ncores = ncores) + computed_pet <- TRUE + } + if (any(is.na(pet[[k]]))) { + mask_na <- which(is.na(pet[[k]])) + warning("There are NAs in PET.") + } + + # Accumulation + diff_p_pet <- data$pr - pet[[k]] + dates_monthly <- .datesmask(dates[[k]]) + accumulated <- .Accumulation(data = diff_p_pet, + dates_monthly = dates_monthly, accum = accum, + time_dim = time_dim, leadtime_dim = leadtime_dim, + ncores = ncores) + + # Standardization + if (standardization) { + spei <- .Standardization(data = accumulated, params = params, + accum = accum, time_dim = time_dim, + leadtime_dim = leadtime_dim, + memb_dim = memb_dim, ref_period = ref_period, + cross_validation = cross_validation, + handle_infinity = handle_infinity, + param_error = param_error, + method = method, distribution = distribution, + na.rm = TRUE, ncores = ncores) + + ref_period <- NULL + params <- spei$params + spei_res <- spei[[1]] + } else { + spei_res <- accumulated + } + pos <- match(names(dim(data[[1]])), names(dim(spei_res))) + spei_res <- aperm(spei_res, pos) + } + + if (standardization) { + if (return_params) { + return(list(spei = spei_res, params = params)) + } else { + return(spei_res) + } + } else { + return(spei_res) + } +} + +.Standardization <- function(data, params = NULL, accum = 1, time_dim = 'syear', + leadtime_dim = 'time', memb_dim = 'ensemble', + ref_period = NULL, cross_validation = FALSE, + handle_infinity = FALSE, param_error = -9999, + method = 'parametric', distribution = 'log-Logistic', + na.rm = FALSE, ncores = NULL) { + + nleadtime <- dim(data)[leadtime_dim] + ntime <- dim(data)[time_dim] + + if (!cross_validation) { + ntime <- 1 + } + + coef = switch(distribution, + "Gamma" = array(NA, dim = 2, dimnames = list(c('alpha', 'beta'))), + "log-Logistic" = array(NA, dim = 3, dimnames = list(c('xi', 'alpha', 'kappa'))), + "PearsonIII" = array(NA, dim = 3, dimnames = list(c('mu', 'sigma', 'gamma')))) + + if (is.null(params)) { + params <- array(NA, dim = c(ntime, nleadtime, length(coef))) + names(dim(params)) <- c(time_dim, leadtime_dim, 'coef') + } else if (length(dim(params)) < 2) { + params <- array(params, dim = c(length(params), ntime, nleadtime)) + # dim(params): [time_dim, leadtime_dim, coef] + # with the values repeated each time_dim and leadtime_dim + params <- aperm(params, c(2,3,1)) + names(dim(params)) <- c(time_dim, leadtime_dim, 'coef') + } + + spei <- Apply(data = list(data = data, params = params), + target_dims = list(data = c(leadtime_dim, time_dim, memb_dim), + params = c(time_dim, leadtime_dim, 'coef')), + fun = .standardization, + coef = coef, + leadtime_dim = leadtime_dim, + time_dim = time_dim, memb_dim = memb_dim, + ref_period = ref_period, handle_infinity = handle_infinity, + cross_validation = cross_validation, param_error = param_error, + method = method, distribution = distribution, + na.rm = na.rm, + output_dims = list(spei = c(leadtime_dim, time_dim, memb_dim), + params = c(time_dim, leadtime_dim, 'coef')), + ncores = ncores) + return(spei) +} + +.standardization <- function(data, params, coef, leadtime_dim = 'time', + time_dim = 'syear', memb_dim = 'ensemble', + ref_period = NULL, handle_infinity = FALSE, + cross_validation = FALSE, param_error = -9999, + method = 'parametric', distribution = 'log-Logistic', + na.rm = FALSE) { + # data: [leadtime_dim, time_dim, memb_dim] + # params: [time_dim, leadtime_dim, 'coef'] + + # maximum number of parameters needed to define any of the considered distributions + ncoef <- length(coef) + nleadtime <- as.numeric(dim(data)[leadtime_dim]) + ntime <- as.numeric(dim(data)[time_dim]) + nmemb <- as.numeric(dim(data)[memb_dim]) + + if (cross_validation) { + params_result <- array(data = NA, dim = c(ntime, nleadtime, ncoef)) + } else { + params_result <- array(data = NA, dim = c(1, nleadtime, ncoef)) + } + names(dim(params_result)) <- c(time_dim, leadtime_dim, 'coef') + + if (all(is.na(data))) { + spei_mod <- array(NA, dim(data)) + # if the data [time, sdate, memb] has no variability it will raise an error + # further down, so we assign a value to the result and skip the step + } else if (anyNA(data) && !na.rm) { + spei_mod <- array(NA, dim(data)) + } else if (var(data, na.rm = T) == 0) { + spei_mod <- array(param_error, dim(data)) # Add this? + } else { + if (is.null(ref_period)) { + ref.start <- NULL + ref.end <- NULL + } else { + ref.start <- ref_period[[1]] + ref.end <- ref_period[[2]] + } + + spei_mod <- array(data = NA, dim = c(nleadtime, ntime, nmemb)) + names(dim(spei_mod)) <- c(leadtime_dim, time_dim, memb_dim) + + for (ff in 1:nleadtime) { + data_subset <- ClimProjDiags::Subset(data, along = leadtime_dim, + indices = ff, drop = 'selected') + params_tmp <- if (all(is.na(params))) {NULL} else {params[, ff, ]} + + spei_data <- .std(data = data_subset, coef = coef, ntime = ntime, + nmemb = nmemb, method = method, + distribution = distribution, na.rm = na.rm, + ref.start = ref.start, ref.end = ref.end, + params = params_tmp, handle_infinity = handle_infinity, + cross_validation = cross_validation) + + spei_mod[ff, , ] <- spei_data[[1]] + params_ff <- spei_data[[2]] + # lengthen dimension coef of params_ff in case it doesn't match the + # corresponding dimension of parms_months + if (!is.null(params_ff)) { + if (length(params_ff) < ncoef) { + params_ff <- append(params_ff, array(NA, dim = ncoef - length(params_ff))) + } + params_result[, ff, ] <- params_ff + } + } + } + return(list(spei = spei_mod, params = params_result)) +} + +.std <- function(data, coef, ntime, nmemb, method = 'parametric', + distribution = 'log-Logistic', na.rm = FALSE, + ref.start = NULL, ref.end = NULL, params = NULL, + handle_infinity = FALSE, cross_validation = FALSE) { + + # data: [time_dim, memb_dim] + # params: NULL or [(ntime), coef] + + fit = 'ub-pwm' # hard-coded + + if (method == 'non-parametric') { + bp <- matrix(0, length(data), 1) + for (i in 1:length(data)) { + bp[i,1] = sum(data[] <= data[i], na.rm = na.rm); # Writes the rank of the data + } + std_index <- qnorm((bp - 0.44)/(length(data) + 0.12)) + dim(std_index) <- c(ntime, nmemb) + # it won't return params to be used in exp_cor; also it is not using + # handle_infinity nor cross_validation + params_result <- NULL + } else { + std_index <- array(NA, c(ntime, nmemb)) + # Select window if necessary + if (!is.null(ref.start) && !is.null(ref.end)) { + data.fit <- window(data, ref.start, ref.end) + } else { + data.fit <- data + } + + if (cross_validation) { + loop_years <- ntime + } else { + loop_years <- 1 + } + + params_result <- array(NA, dim = c(loop_years, length(coef))) + colnames(params_result) <- names(coef) + + for (nsd in 1:loop_years) { # loop over years (in case of cross_validation) + # Cumulative series (acu) + if (cross_validation) { + acu <- as.vector(data.fit[-nsd, ]) + } else { + acu <- as.vector(data.fit) + } + + acu.sorted <- sort.default(acu, method = "quick") + # remove NAs (no need if(na.rm) because if there are NA and na.rm = F + # we don't get to this point) + acu.sorted <- acu.sorted[!is.na(acu.sorted)] + # else all acu was NA and we don't need to continue with this case + + if (length(acu.sorted) != 0) { + acu_sd <- sd(acu.sorted) + if (!is.na(acu_sd)) { + if (acu_sd != 0) { + if (distribution != "log-Logistic") { + pze <- sum(acu == 0) / length(acu) + acu.sorted = acu.sorted[acu.sorted > 0] + } + if (!is.null(params)) { + f_params <- as.vector(params) + params_result[nsd, ] <- f_params + } else { + # else coef will be NA + if (length(acu.sorted) >= 4) { + # Calculate probability weighted moments based on fit with lmomco or TLMoments + pwm = switch(fit, + "pp-pwm" = pwm.pp(acu.sorted, -0.35, 0, nmom = 3), + pwm.ub(acu.sorted, nmom = 3) + # TLMoments::PWM(acu.sorted, order = 0:2) + ) + + # Check L-moments validity + lmom <- pwm2lmom(pwm) + if (!(!are.lmom.valid(lmom) || anyNA(lmom[[1]]) || any(is.nan(lmom[[1]])))) { + # lmom fortran functions need specific inputs L1, L2, T3 + # this is handled by lmomco internally with lmorph + fortran_vec = c(lmom$lambdas[1:2], lmom$ratios[3]) + # Calculate parameters based on distribution with lmom then lmomco + f_params = switch(distribution, + "log-Logistic" = tryCatch(lmom::pelglo(fortran_vec), + error = function(e){parglo(lmom)$para}), + "Gamma" = tryCatch(lmom::pelgam(fortran_vec), + error = function(e){pargam(lmom)$para}), + "PearsonIII" = tryCatch(lmom::pelpe3(fortran_vec), + error = function(e){parpe3(lmom)$para})) + # Adjust if user chose log-Logistic and max-lik + if (distribution == 'log-Logistic' && fit == 'max-lik') { + f_params = parglo.maxlik(acu.sorted, f_params)$para + } + params_result[nsd, ] <- f_params + } # end for the case the L-moments are not valid (std_index will be NA) + } # end for case there are not enough values to estimate the parameters (std_index will be NA) + } # end estimation of f_param + # Calculate cdf based on distribution with lmom + if (all(is.na(params_result[nsd,]))) { + cdf_res <- NA + } else { + f_params <- params_result[nsd,] + f_params <- f_params[which(!is.na(f_params))] + cdf_res = switch(distribution, + "log-Logistic" = lmom::cdfglo(data, f_params), + "Gamma" = lmom::cdfgam(data, f_params), + "PearsonIII" = lmom::cdfpe3(data, f_params)) + } + + std_index_cv <- array(qnorm(cdf_res), dim = c(ntime, nmemb)) + + # Adjust if user chose Gamma or PearsonIII - Not tested: For future development + # if (distribution != 'log-Logistic') { + # std_index[ff,s] = qnorm(pze + (1-pze)*pnorm(std_index[ff,s])) # ff doesn't exist at this point + # } + if (cross_validation) { + std_index[nsd, ] <- std_index_cv[nsd, ] + } else { + std_index <- std_index_cv + } + } + } # end if for the case there is no variability + } # end if for the case all NA in acu + } # next year (in case of cross_validation or all done if cross_validation == F) + + if (handle_infinity) { # could also use "param_error" ?; we are giving it the min/max value of the grid point + std_index[is.infinite(std_index) & std_index < 0] <- min(std_index[!is.infinite(std_index)]) + std_index[is.infinite(std_index) & std_index > 0] <- max(std_index[!is.infinite(std_index)]) + } + # f_params will be params only if cross_validation is FALSE + # (otherwise there will be one f_params per year; + # but the output params will be read only in the case that + # it is called with cross_validation FALSE) + } + return(list(std_index, params_result)) +} + +.Accumulation <- function(data, dates_monthly, accum = 2, + time_dim = 'syear', leadtime_dim = 'time', + ncores = NULL) { + + accumulated <- Apply(data = list(data), + target_dims = list(data = c(leadtime_dim, time_dim)), + dates_monthly = dates_monthly, + accum = accum, + output_dims = c(leadtime_dim, time_dim), + leadtime_dim = leadtime_dim, time_dim = time_dim, + fun = .accumulation, + ncores = ncores)$output1 + + pos <- match(names(dim(accumulated)), names(dim(data))) + accumulated <- aperm(accumulated, pos) + + return(accumulated) + +} + +.accumulation <- function(data, dates_monthly, accum = 1, + time_dim = 'syear', leadtime_dim = 'time') { + + # data:[time, syear] + dims <- dim(data) + + data_vector <- array(NA, dim = length(dates_monthly)) + count <- 1 + for (dd in 1:length(dates_monthly)) { + if (dates_monthly[dd] == 1) { + data_vector[dd] <- as.vector(data)[count] + count <- count + 1 + } + } + + data_sum_x <- rollapply(data_vector, accum, sum) + # adds as many NAs as needed at the begining to account for the months that + # cannot be added (depends on accu) and so that the position in the vector + # corresponds to the accumulated of the previous months (instead of the + # accumulated of the next months) + data_sum_x <- c(rep(NA, accum-1), data_sum_x) + # discard the months that don't appear in the original data + data_sum_x <- data_sum_x[which(dates_monthly == 1)] + + accum_result <- array(data_sum_x, dim = c(dims)) + return(accum_result) +} + +.datesmask <- function(dates) { + ini <- as.Date(paste(lubridate::year(min(dates)), 01, 01, sep = '-')) + end <- as.Date(paste(lubridate::year(max(dates)), 12, 31, sep = '-')) + daily <- as.Date(ini:end) + monthly <- daily[which(lubridate::day(daily) == 1)] + dates_mask <- array(0, dim = length(monthly)) + for (dd in 1:length(dates)) { + ii <- which(monthly == as.Date(paste(lubridate::year(dates[dd]), + lubridate::month(dates[dd]), + 01, sep = '-'))) + dates_mask[ii] <- 1 + } + return(dates_mask) +} \ No newline at end of file diff --git a/R/zzz.R b/R/zzz.R index cf9163970f76c2da74b3d811c1fa0b71445beeab..47d871d4d8d4f2a3d1dfe01b9a79c5aea8674e70 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -88,4 +88,22 @@ wind2CF <- function(wind, pc) { power <- wind2power(wind, pc) CF <- power / pc$attr$RatedPower return(CF) +} + +.KnownLonNames <- function() { + known_lon_names <- c('lon', 'lons', 'longitude', 'x', 'i', 'nav_lon') +} + +.KnownLatNames <- function() { + known_lat_names <- c('lat', 'lats', 'latitude', 'y', 'j', 'nav_lat') +} + +.return2list <- function(data1, data2 = NULL) { + if (is.null(data1) & is.null(data2)) { + return(NULL) + } else if (is.null(data2)) { + return(list(data1)) + } else { + return(list(data1, data2)) + } } \ No newline at end of file diff --git a/man/CST_PeriodPET.Rd b/man/CST_PeriodPET.Rd new file mode 100644 index 0000000000000000000000000000000000000000..10383afcef72d310a54aae7dd30493237d029005 --- /dev/null +++ b/man/CST_PeriodPET.Rd @@ -0,0 +1,85 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/PeriodPET.R +\name{CST_PeriodPET} +\alias{CST_PeriodPET} +\title{Compute the Potential Evapotranspiration} +\usage{ +CST_PeriodPET( + data, + pet_method = "hargreaves", + time_dim = "syear", + leadtime_dim = "time", + lat_dim = "latitude", + na.rm = FALSE, + ncores = NULL +) +} +\arguments{ +\item{data}{A named list with the needed \code{s2dv_cube} objects containing +the seasonal forecast experiment in the data element for each variable. +Specific variables are needed for each method used in computing the +Potential Evapotranspiration. See parameter 'pet_method'. The accepted +variables for each method are: for 'hargreaves', 'tmin' and 'tmax'; for +'hargreaves_modified' are 'tmin', 'tmax' and 'pr'; for method 'thornthwaite' +'tmean' is required. The units for temperature variables +('tmin', 'tmax' and 'tmean') need to be in Celcius degrees; the units for +precipitation ('pr') need to be in mm/month.} + +\item{pet_method}{A character string indicating the method used to compute +the potential evapotranspiration. The accepted methods are: +'hargreaves' and 'hargreaves_modified', that require the data to have +variables tmin and tmax; and 'thornthwaite', that requires variable +'tmean'.} + +\item{time_dim}{A character string indicating the name of the temporal +dimension. By default, it is set to 'syear'.} + +\item{leadtime_dim}{A character string indicating the name of the temporal +dimension. By default, it is set to 'time'.} + +\item{lat_dim}{A character string indicating the name of the latitudinal +dimension. By default it is set by 'latitude'.} + +\item{na.rm}{A logical value indicating whether NA values should be removed +from data. It is FALSE by default.} + +\item{ncores}{An integer value indicating the number of cores to use in +parallel computation.} + +\item{dates}{An array of temporal dimensions containing the Dates of +'data'. It must be of class 'Date' or 'POSIXct'.} + +\item{lat}{A numeric vector containing the latitude values of 'data'.} +} +\description{ +Compute the Potential evapotranspiration (PET) that is the amount of +evaporation and transpiration that would occur if a sufficient water source +were available. This function calculate PET according to the Thornthwaite, +Hargreaves or Hargreaves-modified equations. +} +\details{ +This function is build to work be compatible with other tools in +that work with 's2dv_cube' object class. The input data must be this object +class. If you don't work with 's2dv_cube', see PeriodPET. For more information +on the SPEI calculation, see functions CST_PeriodStandardization and +CST_PeriodAccumulation. +} +\examples{ +dims <- c(time = 3, syear = 3, ensemble = 1, latitude = 1) + +exp_tasmax <- array(rnorm(360, 27.73, 5.26), dim = dims) +exp_tasmin <- array(rnorm(360, 14.83, 3.86), dim = dims) +exp_prlr <- array(rnorm(360, 21.19, 25.64), dim = dims) +end_year <- 2012 +dates_exp <- as.POSIXct(c(paste0(2010:end_year, "-08-16"), + paste0(2010:end_year, "-09-15"), + paste0(2010:end_year, "-10-16")), "UTC") +dim(dates_exp) <- c(syear = 3, time = 3) + +lat <- c(40) + +exp1 <- list('tmax' = exp_tasmax, 'tmin' = exp_tasmin, 'pr' = exp_prlr) + +res <- PeriodPET(data = exp1, lat = lat, dates = dates_exp) + +} diff --git a/man/CST_PeriodSPEI.Rd b/man/CST_PeriodSPEI.Rd new file mode 100644 index 0000000000000000000000000000000000000000..1500a648b3af77766ec6ad85ec3fbad343f24f00 --- /dev/null +++ b/man/CST_PeriodSPEI.Rd @@ -0,0 +1,212 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/PeriodSPEI.R +\name{CST_PeriodSPEI} +\alias{CST_PeriodSPEI} +\title{Compute the Standardised Precipitation-Evapotranspiration Index} +\usage{ +CST_PeriodSPEI( + exp, + exp_cor = NULL, + time_dim = "syear", + leadtime_dim = "time", + memb_dim = "ensemble", + lat_dim = "latitude", + accum = 1, + ref_period = NULL, + params = NULL, + pet_exp = NULL, + pet_expcor = NULL, + standardization = TRUE, + cross_validation = FALSE, + pet_method = "hargreaves", + method = "parametric", + distribution = "log-Logistic", + param_error = -9999, + handle_infinity = FALSE, + return_params = FALSE, + na.rm = FALSE, + ncores = NULL +) +} +\arguments{ +\item{exp}{A named list with the needed \code{s2dv_cube} objects containing +the seasonal forecast experiment in the data element for each variable. +Specific variables are needed for each method used in computing the +Potential Evapotranspiration. See parameter 'pet_method'. The accepted +variables are: 'tmin' and 'tmax', required for methods 'hargreaves' and +'hargreaves_modified' and 'tmean', required for method 'thornthwaite'. +Variable 'pr' is always needed. The units for temperature variables +('tmin', 'tmax' and 'tmean') need to be in Celcius degrees; the units for +precipitation ('pr') need to be in mm/month.} + +\item{exp_cor}{A named list with the needed \code{s2dv_cube} objects for each +variable in which the quantile PeriodSPEI should be applied. If it is not +specified, the PeriodSPEI is calculated from object 'exp'.} + +\item{time_dim}{A character string indicating the name of the temporal +dimension. By default, it is set to 'syear'.} + +\item{leadtime_dim}{A character string indicating the name of the temporal +dimension. By default, it is set to 'time'.} + +\item{memb_dim}{A character string indicating the name of the dimension in +which the ensemble members are stored. When set it to NULL, threshold is +computed for individual members.} + +\item{lat_dim}{A character string indicating the name of the latitudinal +dimension. By default it is set by 'latitude'.} + +\item{accum}{An integer value indicating the number of months for the +accumulation for each variable. When it is greater than 1, the result will +be filled with NA until the accum time_dim dimension number due to the +accumulation to previous months.} + +\item{ref_period}{A list with two numeric values with the starting and end +points of the reference period used for computing the index. The default +value is NULL indicating that the first and end values in data will be +used as starting and end points.} + +\item{params}{A multidimensional array with named dimensions for computing the +SPEI. This option overrides computation of fitting parameters. It needs +to have the following dimensions: same leadtime dimension of 'exp' +(specified in 'leadtime_dim'); time dimension of length 1 (specified in +'time_dim'); and a dimension named 'coef' with the length of the +coefficients needed for the used distribution (for 'Gamma' coef dimension is +of lenght 2, for 'log-Logistic' or 'PearsonIII' is of length 3). It can't +have member dimension (specified in 'memb_dim').} + +\item{pet_exp}{A multidimensional array containing the Potential +EvapoTranspiration data of 'exp'. It must have the same dimensions of the +variable 'pr' of 'exp'. If it is NULL it is calculated using the provided +variables with specified 'pet_method'. It is NULL by default.} + +\item{pet_expcor}{A multidimensional array containing the Potential +EvapoTranspiration data of 'exp_cor'. It must have the same dimensions of +the variable 'pr' of 'exp_cor'. If it is NULL it is calculated using the +provided variables with specified 'pet_method'. It is NULL by default.} + +\item{standardization}{A logical value indicating wether the standardization +is computed.} + +\item{cross_validation}{A logical value indicating if cross validation is +done (TRUE), or not (FALSE). It only will be used when 'exp_cor' and +is not provided. It is FALSE by default.} + +\item{pet_method}{A character string indicating the method used to compute +the potential evapotranspiration. The accepted methods are: +'hargreaves' and 'hargreaves_modified', that require the data to have +variables tmin and tmax; and 'thornthwaite', that requires variable +'tmean'.} + +\item{method}{A character string indicating the standardization method used. +If can be: 'parametric' or 'non-parametric'. It is set to 'parametric' by +default.} + +\item{distribution}{A character string indicating the name of the distribution +function to be used for computing the SPEI. The accepted names are: +'log-Logistic' and 'Gamma'. It is set to 'log-Logistic' by default. The +'Gamma' method only works when only precipitation is provided and other +variables are 0 because it is positive defined (SPI indicator).} + +\item{param_error}{A numeric value with the error accepted.} + +\item{handle_infinity}{A logical value wether to return Infinite values (TRUE) +or not (FALSE).} + +\item{na.rm}{A logical value indicating whether NA values should be removed +from data. It is FALSE by default. If it is FALSE and there are NA values, +(if standardization is TRUE) all values of other dimensions except time_dim +and leadtime_dim will be set to NA directly. On the other hand, if it is +TRUE, if the data from other dimensions except time_dim and leadtime_dim is +not reaching 4 values, it is not enough values to estimate the parameters +and the result will include NA.} + +\item{ncores}{An integer value indicating the number of cores to use in +parallel computation.} +} +\value{ +An 's2dv_cube' object containing the SPEI multidimensional array in +element \code{data}. If 'exp_cor' is provided, only results from 'exp_cor' +will be provided. The parameters of the standardization will only be returned +if 'return_params' is TRUE. The SPEI will only be computed if +'standardization' is TRUE. If 'standardization' is FALSE, only the climatic +water balance (precipitation minus evapotranspiration) will be returned. The +resultant arrays will have the same dimensions as the initial input data. The +other elements in the 's2dv_cube' will be updated with the combined +information of the input data arrays. +} +\description{ +Calculation of the Standardised Precipitation-Evapotranspiration Index (SPEI) +that is a multiscalar drought index based on climatic data. It can be used for +determining the onset, duration and magnitude of drought conditions with +respect to normal conditions in a variety of natural and managed systems such +as crops, ecosystems, rivers, water resources, etc. The idea behind the SPEI +is to compare the highest possible evapotranspiration with the current water +availability. The SPEI uses the monthly (or weekly) difference between +precipitation and potential evapotranspiration. This represents a simple +climatic water balance which is calculated at different time scales to obtain +the SPEI. This function is build to work be compatible with other tools in +that work with 's2dv_cube' object class. The input data must be this object +class. If you don't work with 's2dv_cube', see PeriodSPEI. +} +\details{ +Next, some specifications for the calculation of this indicator will be +discussed. On the one hand, the model to be used to calculate potential +evapotranspiration is specified with the pet_method parameter (hargreaves, +hargraves modified or thornwhite). On the other hand, to choose the time scale +in which you want to accumulate the SPEI (SPEI3, SPEI6...) is done using the +accum parameter, where you must indicate the number of time steps you want to +accumulate throughout leadtime_dim. Since the accumulation is done for the +elapsed time steps, there will be no complete accumulations until reaching the +time instant equal to the value of the parameter. For this reason, in the +result, we will find that for the dimension where the accumulation has been +carried out, the values of the array will be NA since they do not include +complete accumulations. Also, there is a parameter to specify if the +standardization that can be TRUE or FALSE. If it is TRUE, the data is fit to a +probability distribution to transform the original values to standardized +units that are comparable in space and time and at different SPEI time scales. +The na.rm parameter is a logical parameter used to decide whether to remove +the NA values from the data before doing the calculation. It must be taken +into account that if na.rm == FALSE and there is some NA value in the specific +coordinates which the SPEI is computed, standardization cannot be carried out +for those coordinates and therefore, the result will be filled with NA for the +specific coordinates. However, when na.rm == TRUE, if the amount of data for +those specific coordinates is smaller than 4, it will not be possible to carry +out because we will not have enough data and the result will be also filled +with NAs for that coordinates. +} +\examples{ +dims <- c(var = 1, sday = 1, sweek = 1, syear = 6, time = 3, + latitude = 2, longitude = 1, ensemble = 25) +dimscor <- c(var = 1, sday = 1, sweek = 1, syear = 1, time = 3, + latitude = 2, longitude = 1, ensemble = 15) + +dates_exp <- as.POSIXct(c(paste0(2010:2015, "-08-16"), + paste0(2010:2015, "-09-15"), + paste0(2010:2015, "-10-16")), 'UTC') +dim(dates_exp) <- c(sday = 1, sweek = 1, syear = 6, time = 3) +dates_expcor <- as.POSIXct(c(paste0(2020, "-08-16"), paste0(2020, "-09-15"), + paste0(2020, "-10-16")), 'UTC') +dim(dates_expcor) <- c(sday = 1, sweek = 1, syear = 1, time = 3) +lat <- c(40,40.1) + +exp_tmax <- array(rnorm(900, 27.73, 5.26), dim = dims) +exp_tmin <- array(rnorm(900, 14.83, 3.86), dim = dims) +exp_pr <- array(rnorm(900, 21.19, 25.64), dim = dims) + +expcor_tmax <- array(rnorm(540, 29.03, 5.67), dim = dimscor) +expcor_tmin <- array(rnorm(540, 15.70, 4.40), dim = dimscor) +expcor_pr <- array(rnorm(540, 15.62, 21.38), dim = dimscor) + +exp <- list('tmax' = exp_tmax, 'tmin' = exp_tmin, 'pr' = exp_pr) +exp_cor <- list('tmax' = expcor_tmax, 'tmin' = expcor_tmin, + 'pr' = expcor_pr) + +exp <- lapply(exp, CSTools::s2dv_cube, coords = list(latitude = lat), + Dates = dates_exp) +exp_cor <- lapply(exp_cor, CSTools::s2dv_cube, coords = list(latitude = lat), + Dates = dates_expcor) + +res <- CST_PeriodSPEI(exp = exp, exp_cor = exp_cor) + +} diff --git a/man/PeriodPET.Rd b/man/PeriodPET.Rd new file mode 100644 index 0000000000000000000000000000000000000000..7ca1073429f35717fe6a3be22bbf7cffe95e46b0 --- /dev/null +++ b/man/PeriodPET.Rd @@ -0,0 +1,84 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/PeriodPET.R +\name{PeriodPET} +\alias{PeriodPET} +\title{Compute the Potential Evapotranspiration} +\usage{ +PeriodPET( + data, + dates, + lat, + pet_method = "hargreaves", + time_dim = "syear", + leadtime_dim = "time", + lat_dim = "latitude", + na.rm = FALSE, + ncores = NULL +) +} +\arguments{ +\item{data}{A named list with the needed \code{s2dv_cube} objects containing +the seasonal forecast experiment in the data element for each variable. +Specific variables are needed for each method used in computing the +Potential Evapotranspiration. See parameter 'pet_method'. The accepted +variables for each method are: for 'hargreaves', 'tmin' and 'tmax'; for +'hargreaves_modified' are 'tmin', 'tmax' and 'pr'; for method 'thornthwaite' +'tmean' is required. The units for temperature variables +('tmin', 'tmax' and 'tmean') need to be in Celcius degrees; the units for +precipitation ('pr') need to be in mm/month.} + +\item{dates}{An array of temporal dimensions containing the Dates of +'data'. It must be of class 'Date' or 'POSIXct'.} + +\item{lat}{A numeric vector containing the latitude values of 'data'.} + +\item{pet_method}{A character string indicating the method used to compute +the potential evapotranspiration. The accepted methods are: +'hargreaves' and 'hargreaves_modified', that require the data to have +variables tmin and tmax; and 'thornthwaite', that requires variable +'tmean'.} + +\item{time_dim}{A character string indicating the name of the temporal +dimension. By default, it is set to 'syear'.} + +\item{leadtime_dim}{A character string indicating the name of the temporal +dimension. By default, it is set to 'time'.} + +\item{lat_dim}{A character string indicating the name of the latitudinal +dimension. By default it is set by 'latitude'.} + +\item{na.rm}{A logical value indicating whether NA values should be removed +from data. It is FALSE by default.} + +\item{ncores}{An integer value indicating the number of cores to use in +parallel computation.} +} +\description{ +Compute the Potential evapotranspiration (PET) that is the amount of +evaporation and transpiration that would occur if a sufficient water source +were available. This function calculate PET according to the Thornthwaite, +Hargreaves or Hargreaves-modified equations. +} +\details{ +For more information on the SPEI calculation, see functions +PeriodStandardization and PeriodAccumulation. +} +\examples{ +dims <- c(time = 3, syear = 3, ensemble = 1, latitude = 1) + +exp_tasmax <- array(rnorm(360, 27.73, 5.26), dim = dims) +exp_tasmin <- array(rnorm(360, 14.83, 3.86), dim = dims) +exp_prlr <- array(rnorm(360, 21.19, 25.64), dim = dims) +end_year <- 2012 +dates_exp <- as.POSIXct(c(paste0(2010:end_year, "-08-16"), + paste0(2010:end_year, "-09-15"), + paste0(2010:end_year, "-10-16")), "UTC") +dim(dates_exp) <- c(syear = 3, time = 3) + +lat <- c(40) + +exp1 <- list('tmax' = exp_tasmax, 'tmin' = exp_tasmin, 'pr' = exp_prlr) + +res <- PeriodPET(data = exp1, lat = lat, dates = dates_exp) + +} diff --git a/man/PeriodSPEI.Rd b/man/PeriodSPEI.Rd new file mode 100644 index 0000000000000000000000000000000000000000..049e084cd301704f4d157e77b6c96ee64feb4f3e --- /dev/null +++ b/man/PeriodSPEI.Rd @@ -0,0 +1,215 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/PeriodSPEI.R +\name{PeriodSPEI} +\alias{PeriodSPEI} +\title{Compute the Standardised Precipitation-Evapotranspiration Index} +\usage{ +PeriodSPEI( + exp, + dates_exp, + lat, + exp_cor = NULL, + dates_expcor = NULL, + time_dim = "syear", + leadtime_dim = "time", + memb_dim = "ensemble", + lat_dim = "latitude", + accum = 1, + ref_period = NULL, + params = NULL, + pet_exp = NULL, + pet_expcor = NULL, + standardization = TRUE, + cross_validation = FALSE, + pet_method = "hargreaves", + method = "parametric", + distribution = "log-Logistic", + param_error = -9999, + handle_infinity = FALSE, + return_params = FALSE, + na.rm = FALSE, + ncores = NULL +) +} +\arguments{ +\item{exp}{A named list with multidimensional array objects containing +the seasonal forecast experiment in the data element for each variable. +Specific variables are needed for each method used in computing the +Potential Evapotranspiration. See parameter 'pet_method'. The accepted +variables are: 'tmin' and 'tmax', required for methods 'hargreaves' and +'hargreaves_modified' and 'tmean', required for method 'thornthwaite'. +Variable 'pr' is always needed. The units for temperature variables +('tmin', 'tmax' and 'tmean') need to be in Celcius degrees; the units for +precipitation ('pr') need to be in mm/month.} + +\item{dates_exp}{An array of temporal dimensions containing the Dates of +'exp'. It must be of class 'Date' or 'POSIXct'.} + +\item{lat}{A numeric vector containing the latitude values of 'exp'.} + +\item{exp_cor}{A named list with multidimensional array objects for each +variable in which the quantile PeriodSPEI should be applied. If it is not +specified, the PeriodSPEI is calculated from object 'exp'.} + +\item{dates_expcor}{An array of temporal dimensions containing the Dates of +'exp_cor'. It must be of class 'Date' or 'POSIXct'.} + +\item{time_dim}{A character string indicating the name of the temporal +dimension. By default, it is set to 'syear'.} + +\item{leadtime_dim}{A character string indicating the name of the temporal +dimension. By default, it is set to 'time'.} + +\item{memb_dim}{A character string indicating the name of the dimension in +which the ensemble members are stored. When set it to NULL, threshold is +computed for individual members.} + +\item{lat_dim}{A character string indicating the name of the latitudinal +dimension. By default it is set by 'latitude'.} + +\item{accum}{accum An integer value indicating the number of months for the +accumulation for each variable. When it is greater than 1, the result will +be filled with NA until the accum time_dim dimension number due to the +accumulation to previous months.} + +\item{ref_period}{A list with two numeric values with the starting and end +points of the reference period used for computing the index. The default +value is NULL indicating that the first and end values in data will be +used as starting and end points.} + +\item{params}{A multidimensional array with named dimensions for computing the +SPEI. This option overrides computation of fitting parameters. It needs +to be of same time dimension (specified in 'time_dim') of 'exp' and a +dimension named 'coef' with the length of the coefficients needed for the +used distribution (for 'Gamma' coef dimension is of lenght 2, for +'log-Logistic' or 'PearsonIII' is of length 3). It also needs to have a +leadtime dimension (specified in 'leadtime_dim') of length 1.} + +\item{pet_exp}{A multidimensional array containing the Potential +EvapoTranspiration data of 'exp'. It must have the same dimensions of the +variable 'pr' of 'exp'. If it is NULL it is calculated using the provided +variables with specified 'pet_method'. It is NULL by default.} + +\item{pet_expcor}{A multidimensional array containing the Potential +EvapoTranspiration data of 'exp_cor'. It must have the same dimensions of +the variable 'pr' of 'exp_cor'. If it is NULL it is calculated using the +provided variables with specified 'pet_method'. It is NULL by default.} + +\item{standardization}{A logical value indicating wether the standardization +is computed.} + +\item{cross_validation}{A logical value indicating if cross validation is +done (TRUE), or not (FALSE). It only will be used when 'exp_cor' and +is not provided. It is FALSE by default.} + +\item{pet_method}{A character string indicating the method used to compute +the potential evapotranspiration. The accepted methods are: +'hargreaves' and 'hargreaves_modified', that require the data to have +variables tmin and tmax; and 'thornthwaite', that requires variable +'tmean'.} + +\item{method}{A character string indicating the standardization method used. +If can be: 'parametric' or 'non-parametric'.} + +\item{distribution}{A character string indicating the name of the distribution +function to be used for computing the SPEI. The accepted names are: +'log-Logistic' and 'Gamma'. It is set to 'log-Logistic' by default. The +'Gamma' method only works when only precipitation is provided and other +variables are 0 because it is positive defined (SPI indicator).} + +\item{param_error}{A numeric value with the error accepted.} + +\item{handle_infinity}{A logical value wether to return Infinite values (TRUE) +or not (FALSE).} + +\item{return_params}{A logical value indicating wether to return parameters +array (TRUE) or not (FALSE). It is FALSE by default.} + +\item{na.rm}{A logical value indicating whether NA values should be removed +from data. It is FALSE by default. If it is FALSE and there are NA values, +(if standardization is TRUE) all values of other dimensions except time_dim +and leadtime_dim will be set to NA directly. On the other hand, if it is +TRUE, if the data from other dimensions except time_dim and leadtime_dim is +not reaching 4 values, it is not enough values to estimate the parameters +and the result will include NA.} + +\item{ncores}{An integer value indicating the number of cores to use in +parallel computation.} +} +\value{ +An 's2dv_cube' object containing the SPEI multidimensional array in +element \code{data}. If 'exp_cor' is provided, only results from 'exp_cor' +will be provided. The parameters of the standardization will only be returned +if 'return_params' is TRUE. The SPEI will only be computed if +'standardization' is TRUE. If 'standardization' is FALSE, only the climatic +water balance (precipitation minus evapotranspiration) will be returned. The +resultant arrays will have the same dimensions as the initial input data. +} +\description{ +Calculation of the Standardised Precipitation-Evapotranspiration Index (SPEI) +that is a multiscalar drought index based on climatic data. It can be used for +determining the onset, duration and magnitude of drought conditions with +respect to normal conditions in a variety of natural and managed systems such +as crops, ecosystems, rivers, water resources, etc. The idea behind the SPEI +is to compare the highest possible evapotranspiration with the current water +availability. The SPEI uses the monthly (or weekly) difference between +precipitation and potential evapotranspiration. This represents a simple +climatic water balance which is calculated at different time scales to obtain +the SPEI. +} +\details{ +Next, some specifications for the calculation of this indicator will be +discussed. On the one hand, the model to be used to calculate potential +evapotranspiration is specified with the pet_method parameter (hargreaves, +hargraves modified or thornwhite). On the other hand, to choose the time scale +in which you want to accumulate the SPEI (SPEI3, SPEI6...) is done using the +accum parameter, where you must indicate the number of time steps you want to +accumulate throughout leadtime_dim. Since the accumulation is done for the +elapsed time steps, there will be no complete accumulations until reaching the +time instant equal to the value of the parameter. For this reason, in the +result, we will find that for the dimension where the accumulation has been +carried out, the values of the array will be NA since they do not include +complete accumulations. Also, there is a parameter to specify if the +standardization that can be TRUE or FALSE. If it is TRUE, the data is fit to a +probability distribution to transform the original values to standardized +units that are comparable in space and time and at different SPEI time scales. +The na.rm parameter is a logical parameter used to decide whether to remove +the NA values from the data before doing the calculation. It must be taken +into account that if na.rm == FALSE and there is some NA value in the specific +coordinates which the SPEI is computed, standardization cannot be carried out +for those coordinates and therefore, the result will be filled with NA for the +specific coordinates. However, when na.rm == TRUE, if the amount of data for +those specific coordinates is smaller than 4, it will not be possible to carry +out because we will not have enough data and the result will be also filled +with NAs for that coordinates. +} +\examples{ +dims <- c(var = 1, sday = 1, sweek = 1, syear = 6, time = 3, + latitude = 2, longitude = 1, ensemble = 25) +dimscor <- c(var = 1, sday = 1, sweek = 1, syear = 1, time = 3, + latitude = 2, longitude = 1, ensemble = 15) +exp_tmax <- array(rnorm(900, 27.73, 5.26), dim = dims) +exp_tmin <- array(rnorm(900, 14.83, 3.86), dim = dims) +exp_pr <- array(rnorm(900, 21.19, 25.64), dim = dims) + +expcor_tmax <- array(rnorm(540, 29.03, 5.67), dim = dimscor) +expcor_tmin <- array(rnorm(540, 15.70, 4.40), dim = dimscor) +expcor_pr <- array(rnorm(540, 15.62, 21.38), dim = dimscor) + +dates_exp <- as.Date(c(paste0(2010:2015, "-08-16"), + paste0(2010:2015, "-09-15"), + paste0(2010:2015, "-10-16"))) +dim(dates_exp) <- c(sday = 1, sweek = 1, syear = 6, time = 3) +dates_expcor <- as.POSIXct(c(paste0(2020, "-08-16"), + paste0(2020, "-09-15"), + paste0(2020, "-10-16")), 'UTC') +dim(dates_expcor) <- c(sday = 1, sweek = 1, syear = 1, time = 3) +lat <- c(40,40.1) + +exp <- list('tmax' = exp_tmax, 'tmin' = exp_tmin, 'pr' = exp_pr) +exp_cor <- list('tmax' = expcor_tmax, 'tmin' = expcor_tmin, + 'pr' = expcor_pr) +res <- PeriodSPEI(exp = exp, exp_cor = exp_cor, lat = lat, + dates_exp = dates_exp, dates_expcor = dates_expcor) + +} diff --git a/tests/testthat/test-PeriodPET.R b/tests/testthat/test-PeriodPET.R new file mode 100644 index 0000000000000000000000000000000000000000..647d1e08727046185f58b1516d83e2ea920eda6c --- /dev/null +++ b/tests/testthat/test-PeriodPET.R @@ -0,0 +1,197 @@ +############################################## + +# cube1 +cube1 <- NULL +cube1$data <- array(rnorm(10), dim = c(lat = 2, time = 5)) +class(cube1) <- 's2dv_cube' + +# cube2 +cube2 <- NULL +cube2$data <- array(rnorm(10), dim = c(lat = 2, time = 5)) +class(cube2) <- 's2dv_cube' +cube2$coords <- list(lat = 1:2) + +# dat1 +dims <- c(syear = 6, time = 3, latitude = 2, longitude = 1, ensemble = 10) + +set.seed(1) +exp_tmax <- array(rnorm(360, 27.73, 5.26), dim = dims) +set.seed(2) +exp_tmin <- array(rnorm(360, 14.83, 3.86), dim = dims) +set.seed(3) +exp_pr <- array(rnorm(360, 21.19, 25.64), dim = dims) + +dates_exp <- as.Date(c(paste0(2010:2015, "-08-16"), paste0(2010:2015, "-09-15"), + paste0(2010:2015, "-10-16"))) +dim(dates_exp) <- c(syear = 6, time = 3) + +lat <- c(40,40.1) + +exp1 <- list('tmax' = exp_tmax, 'tmin' = exp_tmin, 'pr' = exp_pr) + +# dat2 +dims2 <- c(styear = 6, ftime = 3, lat = 2, lon = 1, member = 10) + +set.seed(1) +exp_tmean <- array(rnorm(100, 17.34, 9.18), dim = dims2) +set.seed(2) +exp_pr <- array(rnorm(360, 21.19, 25.64), dim = dims2) + +dates_exp2 <- as.Date(c(paste0(2010:2015, "-08-16"), paste0(2010:2015, "-09-15"), + paste0(2010:2015, "-10-16"))) +dim(dates_exp2) <- c(sday = 1, sweek = 1, styear = 6, ftime = 3) + +lat <- c(40,40.1) + +exp2 <- list('tmean' = exp_tmean, 'pr' = exp_pr) + +# cube4 +cube4_exp <- lapply(exp1, function(x) { + suppressWarnings( + CSTools::s2dv_cube(data = x, coords = list(latitude = c(40, 40.1)), + varName = 'test', Dates = as.POSIXct(dates_exp)) + ) +}) + +############################################## + +test_that("1. Initial checks CST_PeriodPET", { + # Check 's2dv_cube' + expect_error( + CST_PeriodPET(data = NULL), + "Parameter 'data' cannot be NULL." + ) + expect_error( + CST_PeriodPET(data = array(10)), + "Parameter 'data' must be a list of 's2dv_cube' class." + ) + # latitude + expect_error( + CST_PeriodPET(data = list(cube1)), + paste0("Spatial coordinate names of parameter 'data' do not match any ", + "of the names accepted by the package.") + ) + # Dates + expect_error( + CST_PeriodPET(data = list(cube2)), + paste0("Element 'Dates' is not found in 'attrs' list of 'data'. ", + "See 's2dv_cube' object description in README file for more ", + "information.") + ) +}) + +############################################## + +test_that("1. Initial checks PeriodPET", { + # data + expect_error( + PeriodPET(data = NULL), + "Parameter 'data' needs to be a named list with the needed variables." + ) + expect_error( + PeriodPET(data = list(1)), + "Parameter 'data' needs to be a named list with the variable names." + ) + expect_error( + PeriodPET(data = list(tmax = array(10))), + "Parameter 'data' needs to be a list of arrays with dimension names." + ) + expect_error( + PeriodPET(data = list(tmax = array(10, c(time = 10)), + tmin = array(10, c(time = 11)))), + "Parameter 'data' variables need to have the same dimensions." + ) + expect_error( + PeriodPET(data = list(tmax = array(10, c(time = 10)), + tmin = array(10, c(ftime = 10)))), + "Parameter 'data' variables need to have the same dimensions." + ) + # lat + expect_error( + PeriodPET(data = exp1, lat = 'lat'), + "Parameter 'lat' must be numeric." + ) + expect_error( + PeriodPET(data = list(tmax = array(10, c(time = 10)), + tmin = array(10, c(time = 10))), lat = 1:2), + "Parameter 'data' must have 'lat_dim' dimension." + ) + # data (2) + expect_warning( + PeriodPET(data = exp1, pet_method = '1', dates = dates_exp, lat = lat), + paste0("Parameter 'pet_method' needs to be 'hargreaves' or ", + "'hargreaves_modified'. It is set to 'hargreaves_modified'.") + ) + # time_dim + expect_error( + PeriodPET(data = exp1, time_dim = 1, dates = dates_exp, lat = lat), + "Parameter 'time_dim' must be a character string." + ) + expect_error( + PeriodPET(data = exp2, lat = lat, dates = dates_exp2, + lat_dim = 'lat', pet_method = 'thornthwaite'), + "Parameter 'time_dim' is not found in 'data' dimension." + ) + # leadtime_dim + expect_error( + PeriodPET(data = exp1, leadtime_dim = 1, dates = dates_exp, lat = lat), + "Parameter 'leadtime_dim' must be a character string." + ) + expect_error( + PeriodPET(data = exp2, lat = lat, dates = dates_exp2, + lat_dim = 'lat', time_dim = 'ftime', pet_method = 'thornthwaite'), + "Parameter 'leadtime_dim' is not found in 'data' dimension." + ) + # lat_dim + expect_error( + PeriodPET(data = exp1, lat_dim = 1, dates = dates_exp, lat = lat) + ) + expect_error( + PeriodPET(data = exp2, lat = lat, dates = dates_exp2), + "Parameter 'data' must have 'lat_dim' dimension." + ) + # na.rm + expect_error( + PeriodPET(data = exp1, na.rm = 1.5, dates = dates_exp, lat = lat), + "Parameter 'na.rm' must be one logical value." + ) + # ncores + expect_error( + PeriodPET(data = exp1, ncores = 1.5, dates = dates_exp, lat = lat), + "Parameter 'ncores' must be a positive integer." + ) +}) + +############################################## + +test_that("2. Output checks", { + res01 <- CST_PeriodPET(data = cube4_exp) + res1 <- PeriodPET(data = exp1, lat = lat, dates = dates_exp) + res2 <- PeriodPET(data = exp2, lat = lat, dates = dates_exp2, + pet_method = c('thornthwaite'), + lat_dim = 'lat', time_dim = 'styear', + leadtime_dim = 'ftime') + # structure + expect_equal( + names(res01), + c('data', 'dims', 'coords', 'attrs') + ) + # dims + expect_equal( + dim(res1), + c(syear = 6, time = 3, latitude = 2, longitude = 1, ensemble = 10) + ) + # values + expect_equal( + res1[1:4], + c(137.77342, 154.55548, 65.72859, 222.20438), + tolerance = 0.0001 + ) + expect_equal( + res2[1:4], + c(77.76124, 118.94212, 66.57568, 185.67074), + tolerance = 0.0001 + ) +}) + +############################################## \ No newline at end of file diff --git a/tests/testthat/test-PeriodSPEI.R b/tests/testthat/test-PeriodSPEI.R new file mode 100644 index 0000000000000000000000000000000000000000..2f1d7b0745bf4e507fb094965f6abf50ff4b7e10 --- /dev/null +++ b/tests/testthat/test-PeriodSPEI.R @@ -0,0 +1,481 @@ +############################################## +# cube1 +cube1 <- NULL +cube1$data <- array(rnorm(10), dim = c(lat = 2, time = 5)) +class(cube1) <- 's2dv_cube' + +# cube2 +cube2 <- NULL +cube2$data <- array(rnorm(10), dim = c(lat = 2, time = 5)) +class(cube2) <- 's2dv_cube' +cube2$coords <- list(lat = 1:2) + +# cube3 +cube3 <- NULL +cube3$data <- array(rnorm(10), dim = c(lat = 2, time = 5)) +class(cube3) <- 's2dv_cube' +cube3$coords <- list(lat = 1:2) +cube3$attrs$Dates <- as.Date(c(paste0(2010:2014, "-08-16"))) + +# dat1 +dims <- c(syear = 6, time = 3, latitude = 2, longitude = 1, ensemble = 10) +dimscor <- c(syear = 1, time = 3, latitude = 2, longitude = 1, ensemble = 15) + +set.seed(1) +exp_tmax <- array(rnorm(360, 27.73, 5.26), dim = dims) +set.seed(2) +exp_tmin <- array(rnorm(360, 14.83, 3.86), dim = dims) +set.seed(3) +exp_pr <- array(rnorm(360, 21.19, 25.64), dim = dims) + +set.seed(1) +expcor_tmax <- array(rnorm(60, 29.03, 5.67), dim = dimscor) +set.seed(2) +expcor_tmin <- array(rnorm(60, 15.70, 4.40), dim = dimscor) +set.seed(3) +expcor_pr <- array(rnorm(60, 15.62, 21.38), dim = dimscor) + +dates_exp <- as.POSIXct(c(paste0(2010:2015, "-08-16"), paste0(2010:2015, "-09-15"), + paste0(2010:2015, "-10-16")), "UTC") +dim(dates_exp) <- c(syear = 6, time = 3) + +dates_expcor <- as.POSIXct(c(paste0(2020, "-08-16"), paste0(2020, "-09-15"), + paste0(2020, "-10-16")), "UTC") +dim(dates_expcor) <- c(syear = 1, time = 3) + +lat <- c(40,40.1) + +exp1 <- list('tmax' = exp_tmax, 'tmin' = exp_tmin, 'pr' = exp_pr) +exp_cor1 <- list('tmax' = expcor_tmax, 'tmin' = expcor_tmin, 'pr' = expcor_pr) +params1 <- array(abs(rnorm(100)), dim = c(syear = 1, time = 3, latitude = 2, + longitude = 1, coef = 3)) + +# dat2 +dims2 <- c(styear = 6, ftime = 3, lat = 2, lon = 1, member = 10) +dimscor2 <- c(styear = 1, ftime = 3, lat = 2, lon = 1, member = 15) + +set.seed(1) +exp_tmean <- array(rnorm(100, 17.34, 9.18), dim = dims2) +set.seed(2) +exp_pr <- array(rnorm(360, 21.19, 25.64), dim = dims2) + +set.seed(1) +expcor_tmean <- array(rnorm(100, 17.23, 9.19), dim = dimscor2) +set.seed(2) +expcor_pr <- array(rnorm(60, 15.62, 21.38), dim = dimscor2) + +dates_exp2 <- as.Date(c(paste0(2010:2015, "-08-16"), paste0(2010:2015, "-09-15"), + paste0(2010:2015, "-10-16"))) +dim(dates_exp2) <- c(sday = 1, sweek = 1, styear = 6, ftime = 3) + +dates_expcor2 <- as.Date(c(paste0(2020, "-08-16"), paste0(2020, "-09-15"), + paste0(2020, "-10-16"))) +dim(dates_expcor2) <- c(sday = 1, sweek = 1, styear = 1, ftime = 3) + +lat <- c(40,40.1) + +exp2 <- list('tmean' = exp_tmean, 'pr' = exp_pr) +exp_cor2 <- list('tmean' = expcor_tmean, 'pr' = expcor_pr) + +# cube4 +cube4_exp <- lapply(exp1, function(x) { + suppressWarnings( + CSTools::s2dv_cube(data = x, coords = list(latitude = c(40, 40.1)), + varName = 'test', Dates = dates_exp) + ) +}) +cube4_expcor <- lapply(exp_cor1, function(x) { + suppressWarnings( + CSTools::s2dv_cube(data = x, coords = list(latitude = c(40, 40.1)), + varName = 'test', Dates = dates_expcor) + ) +}) + +############################################## + +test_that("1. Initial checks CST_PeriodSPEI", { + # Check 's2dv_cube' + expect_error( + CST_PeriodSPEI(exp = NULL), + "Parameter 'exp' cannot be NULL." + ) + expect_error( + CST_PeriodSPEI(exp = array(10)), + "Parameter 'exp' must be a list of 's2dv_cube' class." + ) + # latitude + expect_error( + CST_PeriodSPEI(exp = list(cube1)), + paste0("Spatial coordinate names of parameter 'exp' do not match any ", + "of the names accepted by the package.") + ) + # Dates + expect_error( + CST_PeriodSPEI(exp = list(cube2)), + paste0("Element 'Dates' is not found in 'attrs' list of 'exp'. ", + "See 's2dv_cube' object description in README file for more ", + "information.") + ) + expect_error( + CST_PeriodSPEI(exp = list(cube3), exp_cor = list(cube2)), + paste0("Element 'Dates' is not found in 'attrs' list of 'exp_cor'. ", + "See 's2dv_cube' object description in README file for more ", + "information.") + ) +}) + +############################################## + +test_that("1. Initial checks PeriodSPEI", { + # exp + expect_error( + PeriodSPEI(exp = NULL), + "Parameter 'exp' needs to be a named list with the needed variables." + ) + expect_error( + PeriodSPEI(exp = list(1)), + "Parameter 'exp' needs to be a named list with the variable names." + ) + expect_error( + PeriodSPEI(exp = list(tmax = array(10))), + "Parameter 'exp' needs to be a list of arrays with dimension names." + ) + expect_error( + PeriodSPEI(exp = list(tmax = array(10, c(time = 10)), + tmin = array(10, c(time = 11)))), + "Parameter 'exp' variables need to have the same dimensions." + ) + expect_error( + PeriodSPEI(exp = list(tmax = array(10, c(time = 10)), + tmin = array(10, c(ftime = 10)))), + "Parameter 'exp' variables need to have the same dimensions." + ) + # exp_cor + expect_error( + PeriodSPEI(exp = exp1, exp_cor = 1), + paste0("Parameter 'exp_cor' needs to be a named list with the needed ", + "variables if it is not NULL.") + ) + expect_error( + PeriodSPEI(exp = exp1, exp_cor = list(1)), + "Parameter 'exp_cor' needs to be a named list with the variable names." + ) + expect_error( + PeriodSPEI(exp = exp1, exp_cor = list('tmax' = array(10))), + "Parameter 'exp_cor' needs to be a list of arrays with dimension names." + ) + expect_error( + PeriodSPEI(exp = exp1, exp_cor = list(tmax = array(10, c(time = 10)), + tmin = array(10, c(time = 11)))), + "Parameter 'exp_cor' variables need to have the same dimensions." + ) + expect_error( + PeriodSPEI(exp = exp1, lat = 'lat', dates_exp = dates_exp), + "Parameter 'lat' must be numeric." + ) + expect_error( + PeriodSPEI(exp = list(pr = array(10, c(time = 10, syear = 1, ensemble = 1))), + lat = 1:2, dates_exp = dates_exp), + "Parameter 'lat_dim' is not found in 'exp' dimension." + ) + # exp (2) + expect_warning( + PeriodSPEI(exp = exp1, pet_method = '1', dates_exp = dates_exp, lat = lat), + paste0("Parameter 'pet_method' needs to be 'hargreaves' or ", + "'hargreaves_modified'. It is set to 'hargreaves_modified'.") + ) + # time_dim + expect_error( + PeriodSPEI(exp = exp1, time_dim = 1, dates_exp = dates_exp, lat = lat), + "Parameter 'time_dim' must be a character string." + ) + expect_error( + PeriodSPEI(exp = exp2, exp_cor = exp_cor2, lat = lat, + dates_exp = dates_exp2, dates_expcor = dates_expcor2, + lat_dim = 'lat', pet_method = 'thornthwaite'), + "Parameter 'time_dim' is not found in 'exp' dimension." + ) + # leadtime_dim + expect_error( + PeriodSPEI(exp = exp1, leadtime_dim = 1, dates_exp = dates_exp, lat = lat), + "Parameter 'leadtime_dim' must be a character string." + ) + expect_error( + PeriodSPEI(exp = exp2, exp_cor = exp_cor2, lat = lat, + dates_exp = dates_exp2, dates_expcor = dates_expcor2, + lat_dim = 'lat', time_dim = 'ftime', pet_method = 'thornthwaite'), + "Parameter 'leadtime_dim' is not found in 'exp' dimension." + ) + # memb_dim + expect_error( + PeriodSPEI(exp = exp1, memb_dim = 1, dates_exp = dates_exp, lat = lat), + "Parameter 'memb_dim' must be a character string." + ) + expect_error( + PeriodSPEI(exp = exp2, exp_cor = exp_cor2, lat = lat, + dates_exp = dates_exp2, dates_expcor = dates_expcor2, + lat_dim = 'lat', time_dim = 'ftime', leadtime_dim = 'styear', + pet_method = 'thornthwaite'), + "Parameter 'memb_dim' is not found in 'exp' dimension." + ) + # lat_dim + expect_error( + PeriodSPEI(exp = exp1, lat_dim = 1, dates_exp = dates_exp, lat = lat) + ) + expect_error( + PeriodSPEI(exp = exp2, exp_cor = exp_cor2, lat = lat, + dates_exp = dates_exp2, dates_expcor = dates_expcor2), + "Parameter 'time_dim' is not found in 'exp' dimension." + ) + # accum + expect_error( + PeriodSPEI(exp = exp1, accum = 10, dates_exp = dates_exp, lat = lat), + "Cannot compute accumulation of 10 months because loaded data has only 3 months." + ) + # ref_period + expect_warning( + PeriodSPEI(exp = exp1, exp_cor = exp_cor1, dates_exp = dates_exp, + dates_expcor = dates_expcor, lat = lat, ref_period = 1), + paste0("Parameter 'ref_period' must be of length two indicating the ", + "first and end years of the reference period. It will not ", + "be used.") + ) + expect_warning( + PeriodSPEI(exp = exp1, exp_cor = exp_cor1, dates_exp = dates_exp, + dates_expcor = dates_expcor, lat = lat, ref_period = list('a', 1)), + paste0("Parameter 'ref_period' must be a numeric vector indicating the ", + "'start' and 'end' years of the reference period. It will not ", + "be used.") + ) + expect_warning( + PeriodSPEI(exp = exp1, exp_cor = exp_cor1, dates_exp = dates_exp, + dates_expcor = dates_expcor, lat = lat, ref_period = list(2012, 2011)), + paste0("In parameter 'ref_period' 'start' cannot be after 'end'. It ", + "will not be used.") + ) + expect_warning( + PeriodSPEI(exp = exp1, exp_cor = exp_cor1, dates_exp = dates_exp, + dates_expcor = dates_expcor, lat = lat, ref_period = list(2008, 2021)), + paste0("Parameter 'ref_period' contain years outside the dates. ", + "It will not be used.") + ) + # standardization + expect_error( + PeriodSPEI(exp = exp1, standardization = 10, dates_exp = dates_exp, lat = lat), + "Parameter 'standardization' must be a logical value." + ) + # param_error + expect_error( + PeriodSPEI(exp = exp1, param_error = TRUE, dates_exp = dates_exp, lat = lat), + "Parameter 'param_error' must be a numeric value." + ) + # handle_infinity + expect_error( + PeriodSPEI(exp = exp1, handle_infinity = 1, dates_exp = dates_exp, lat = lat), + "Parameter 'handle_infinity' must be a logical value." + ) + # cross_validation + expect_error( + PeriodSPEI(exp = exp1, cross_validation = 1, dates_exp = dates_exp, lat = lat), + "Parameter 'cross_validation' must be a logical value." + ) + # method + expect_error( + PeriodSPEI(exp = exp1, method = 1, dates_exp = dates_exp, lat = lat), + paste0("Parameter 'method' must be a character string containing one ", + "of the following methods: 'parametric' or 'non-parametric'.") + ) + # distribution + expect_error( + PeriodSPEI(exp = exp1, distribution = 1, dates_exp = dates_exp, lat = lat), + paste0("Parameter 'distribution' must be a character string containing one ", + "of the following distributions: 'log-Logistic', 'Gamma' or ", + "'PearsonIII'.") + ) + # ncores + expect_error( + PeriodSPEI(exp = exp1, ncores = 1.5, dates_exp = dates_exp, lat = lat), + "Parameter 'ncores' must be a positive integer." + ) +}) + +############################################## + +test_that("2. Output checks: CST_PeriodSPEI", { + res1 <- CST_PeriodSPEI(exp = cube4_exp, exp_cor = NULL) + res2 <- CST_PeriodSPEI(exp = cube4_exp, exp_cor = cube4_expcor) + res3 <- CST_PeriodSPEI(exp = cube4_exp, exp_cor = cube4_expcor, standardization = F) + res4 <- CST_PeriodSPEI(exp = cube4_exp, exp_cor = NULL, return_params = T) + expect_equal( + names(res1), + c("data", "dims", "coords", "attrs") + ) + expect_equal( + res2$attrs$Variable$varName, + "SPEI" + ) + expect_equal( + res3$attrs$Variable$varName, + "Precipitation minus accumulated PET" + ) + expect_equal( + names(res3), + c("data", "dims", "coords", "attrs") + ) + expect_equal( + names(res4), + c("spei", "params") + ) +}) + +############################################## + +test_that("2. Output checks", { + res1 <- PeriodSPEI(exp = exp1, exp_cor = exp_cor1, lat = lat, + dates_exp = dates_exp, dates_expcor = dates_expcor, + return_params = TRUE) + res2 <- PeriodSPEI(exp = exp1, exp_cor = exp_cor1, lat = lat, + dates_exp = dates_exp, dates_expcor = dates_expcor, + standardization = FALSE) # No info about accumulation + res3 <- PeriodSPEI(exp = exp1, exp_cor = NULL, lat = lat, + dates_exp = dates_exp, return_params = TRUE) + # output dims + expect_equal( + names(res1), + c('spei', 'params') + ) + expect_equal( + dim(res2), + c(syear = 1, time = 3, latitude = 2, longitude = 1, ensemble = 15) + ) + expect_equal( + names(res3), + c('spei', 'params') + ) + expect_equal( + dim(res1[[1]]), + dimscor + ) + expect_equal( + dim(res1[[2]])[which(!names(dim(res1[[2]])) %in% c('coef', 'syear'))], + dims[which(!names(dims) %in% c('syear', 'ensemble'))] + ) + expect_equal( + dim(res2), + dimscor + ) + expect_equal( + dim(res3[[2]]), + c(syear = 1, time = 3, coef = 3, latitude = 2, longitude = 1) + ) + # exp + # exp_cor + # pet + # accum + res11 <- PeriodSPEI(exp = exp1, exp_cor = NULL, lat = lat, accum = 2, + dates_exp = dates_exp, na.rm = TRUE) + # expect_equal( + # res11[1,3,1,1,][1:4], + # c(-0.4292409, -0.1375149, -0.5564081, -0.4273380), + # tolerance = 0.0001 + # ) + # ref_period + res_ref <- PeriodSPEI(exp = exp1, exp_cor = exp_cor1, lat = lat, accum = 2, + dates_exp = dates_exp, dates_expcor = dates_expcor, + na.rm = TRUE, ref_period = list(2011, 2013)) + expect_equal( + !identical(res1[[1]], res_ref), + TRUE + ) + # params + res5 <- PeriodSPEI(exp = exp1, lat = lat, dates_exp = dates_exp, + params = params1, return_params = TRUE) + expect_error( + PeriodSPEI(exp = exp1, lat = lat, dates_exp = dates_exp, + params = array(abs(rnorm(100)), dim = dimscor)), + paste0("Parameter 'params' must be a multidimensional array with named ", + "dimensions: 'syear', 'time' and 'coef'.") + ) + expect_equal( + dim(res5$params), + c(syear = 1, time = 3, coef = 3, latitude = 2, longitude = 1) + ) + # standarization + res4 <- PeriodSPEI(exp = exp1, exp_cor = NULL, lat = lat, + dates_exp = dates_exp, standardization = FALSE) + expect_equal( + names(res4), + NULL + ) + expect_equal( + dim(res4), + c(syear = 6, time = 3, latitude = 2, longitude = 1, ensemble = 10) + ) + # cross_validation + expect_warning( + PeriodSPEI(exp = exp1, exp_cor = exp_cor1, lat = lat, + dates_exp = dates_exp, dates_expcor = dates_expcor, + cross_validation = TRUE), + paste0("Detected 'cross_validation' = TRUE. It will be set as FALSE ", + "since 'exp_cor' is provided.") + ) + res_crossval_T <- PeriodSPEI(exp = exp1, lat = lat, dates_exp = dates_exp, + cross_validation = TRUE, return_params = TRUE) + res_crossval_F <- PeriodSPEI(exp = exp1, lat = lat, dates_exp = dates_exp, + cross_validation = FALSE, return_params = TRUE) + # cross_validation = TRUE + expect_equal( + dim(res_crossval_T$spei), + dims + ) + expect_equal( + dim(res_crossval_T$params), + c(syear = 6, time = 3, coef = 3, latitude = 2, longitude = 1) + ) + # cross_validation = FALSE + expect_equal( + dim(res_crossval_F$params)[-which(names(dim(res_crossval_F$params)) == 'coef')], + dimscor[-which(names(dimscor) == 'ensemble')] + ) + + # pet_method - ok + res5 <- PeriodSPEI(exp = exp1, exp_cor = exp_cor1, lat = lat, + dates_exp = dates_exp, dates_expcor = dates_expcor, + pet_method = c('hargreaves', 'hargreaves_modified')) + res6 <- PeriodSPEI(exp = exp1, exp_cor = exp_cor1, lat = lat, + dates_exp = dates_exp, dates_expcor = dates_expcor, + pet_method = c('hargreaves_modified', 'hargreaves')) + expect_equal( + identical(res5, res6), + FALSE + ) + expect_equal( + dim(res5), + dim(res6) + ) + + # time_dim, leadtime_dim, memb_dim, lat_dim + res7 <- PeriodSPEI(exp = exp2, exp_cor = exp_cor2, lat = lat, + dates_exp = dates_exp2, dates_expcor = dates_expcor2, + pet_method = c('thornthwaite', 'thornthwaite'), + lat_dim = 'lat', time_dim = 'styear', + leadtime_dim = 'ftime', memb_dim = 'member') + # method - ok + res8 <- PeriodSPEI(exp = exp1, exp_cor = exp_cor1, lat = lat, + dates_exp = dates_exp, dates_expcor = dates_expcor, + method = 'non-parametric') + # distribution - Only works for 'log-Logistic' + res9 <- PeriodSPEI(exp = exp1, exp_cor = NULL, lat = lat, + dates_exp = dates_exp, distribution = 'PearsonIII') # NA + res10 <- PeriodSPEI(exp = exp1, exp_cor = NULL, lat = lat, + dates_exp = dates_exp, distribution = 'Gamma') # NA + # param_error - + # handle_infinity - OK + res9 <- PeriodSPEI(exp = exp1, exp_cor = NULL, lat = lat, + dates_exp = dates_exp, handle_infinity = FALSE) + # na.rm - + + # ncores +}) + +############################################## \ No newline at end of file