From 5a4f6763d4a8284cc1e91b2f5fddd01d367d660b Mon Sep 17 00:00:00 2001 From: aho Date: Mon, 2 Nov 2020 12:02:54 +0100 Subject: [PATCH 1/3] Update Load() and Utils.R from s2dverification --- R/Load.R | 198 +++++++++++++++++++++++++++++++++++++++++++----------- R/Utils.R | 61 +++++------------ 2 files changed, 175 insertions(+), 84 deletions(-) diff --git a/R/Load.R b/R/Load.R index b6a3ef6..9feec5d 100644 --- a/R/Load.R +++ b/R/Load.R @@ -42,7 +42,7 @@ #'the individual grid of each dataset but can also be averaged after #'interpolating into a common grid. See parameters 'grid' and 'method'.\cr #'Once the two arrays are filled by calling this function, other functions in -#'the s2dv package that receive as inputs data formatted in this +#'the s2dverification package that receive as inputs data formatted in this #'data structure can be executed (e.g: \code{Clim()} to compute climatologies, #'\code{Ano()} to compute anomalies, ...).\cr\cr #'Load() has many additional parameters to disable values and trim dimensions @@ -455,7 +455,7 @@ #' Warning: list() compulsory even if loading 1 experimental dataset only!\cr #' E.g., list(array(1, dim = c(num_lons, num_lats))) #'@param maskobs See help on parameter 'maskmod'. -#'@param configfile Path to the s2dv configuration file from which +#'@param configfile Path to the s2dverification configuration file from which #' to retrieve information on location in file system (and other) of datasets.\cr #' If not specified, the configuration file used at BSC-ES will be used #' (it is included in the package).\cr @@ -670,6 +670,13 @@ #' } #' } #' +#'@keywords datagen +#'@author History:\cr +#'0.1 - 2011-03 (V. Guemas, \email{virginie.guemas@bsc.es}) - Original code\cr +#'1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens@bsc.es}) - Formatting to CRAN\cr +#'1.2 - 2015-02 (N. Manubens, \email{nicolau.manubens@bsc.es}) - Generalisation + parallelisation\cr +#'1.3 - 2015-07 (N. Manubens, \email{nicolau.manubens@bsc.es}) - Improvements related to configuration file mechanism\cr +#'1.4 - 2016-01 (N. Manubens, \email{nicolau.manubens@bsc.es}) - Added subsetting capabilities #'@examples #'# Let's assume we want to perform verification with data of a variable #'# called 'tos' from a model called 'model' and observed data coming from @@ -755,7 +762,7 @@ #'# Example 1: Providing lists of lists to 'exp' and 'obs': #'# #' \dontrun{ -#'data_path <- system.file('sample_data', package = 's2dv') +#'data_path <- system.file('sample_data', package = 's2dverification') #'exp <- list( #' name = 'experiment', #' path = file.path(data_path, 'model/$EXP_NAME$/monthly_mean', @@ -781,7 +788,7 @@ #'# writing a configuration file). #'# #' \dontrun{ -#'data_path <- system.file('sample_data', package = 's2dv') +#'data_path <- system.file('sample_data', package = 's2dverification') #'expA <- list(name = 'experiment', path = file.path(data_path, #' 'model/$EXP_NAME$/$STORE_FREQ$_mean/$VAR_NAME$_3hourly', #' '$VAR_NAME$_$START_DATE$.nc')) @@ -795,7 +802,7 @@ #' output = 'areave', latmin = 27, latmax = 48, #' lonmin = -12, lonmax = 40) #'# -#'# Example 3: providing character strings in 'exp' and 'obs', and providing +#'# Example 2: providing character strings in 'exp' and 'obs', and providing #'# a configuration file. #'# The configuration file 'sample.conf' that we will create in the example #'# has the proper entries to load these (see ?LoadConfigFile for details on @@ -806,7 +813,7 @@ #'c <- ConfigFileOpen(configfile) #'c <- ConfigEditDefinition(c, 'DEFAULT_VAR_MIN', '-1e19', confirm = FALSE) #'c <- ConfigEditDefinition(c, 'DEFAULT_VAR_MAX', '1e19', confirm = FALSE) -#'data_path <- system.file('sample_data', package = 's2dv') +#'data_path <- system.file('sample_data', package = 's2dverification') #'exp_data_path <- paste0(data_path, '/model/$EXP_NAME$/') #'obs_data_path <- paste0(data_path, '/$OBS_NAME$/') #'c <- ConfigAddEntry(c, 'experiments', dataset_name = 'experiment', @@ -825,14 +832,15 @@ #' } #' \dontshow{ #'startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101') -#'sampleData <- s2dv:::.LoadSampleData('tos', c('experiment'), -#' c('observation'), startDates, -#' output = 'areave', -#' latmin = 27, latmax = 48, -#' lonmin = -12, lonmax = 40) +#'sampleData <- s2dverification:::.LoadSampleData('tos', c('experiment'), +#' c('observation'), startDates, +#' output = 'areave', +#' latmin = 27, latmax = 48, +#' lonmin = -12, lonmax = 40) #' } #'@import parallel bigmemory methods #'@importFrom stats ts window na.omit +#'@importFrom abind abind #'@export Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, nmemberobs = NULL, nleadtime = NULL, leadtimemin = 1, @@ -843,6 +851,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, configfile = NULL, varmin = NULL, varmax = NULL, silent = FALSE, nprocs = NULL, dimnames = NULL, remapcells = 2, path_glob_permissive = 'partial') { + #library(parallel) #library(bigmemory) @@ -1045,7 +1054,6 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, if (!is.character(storefreq) || !(storefreq %in% c('monthly', 'daily'))) { stop("Error: parameter 'storefreq' is wrong, can take value 'daily' or 'monthly'.") } - # sampleperiod if (is.null(sampleperiod) || !is.numeric(sampleperiod)) { stop("Error: parameter 'sampleperiod' is wrong. It should be numeric.") @@ -1163,9 +1171,9 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, # configfile if (is.null(configfile)) { - configfile <- system.file("config", "BSC.conf", package = "s2dv") + configfile <- system.file("config", "BSC.conf", package = "s2dverification") } else if (!is.character(configfile) || !(nchar(configfile) > 0)) { - stop("Error: parameter 'configfile' must be a character string with the path to an s2dv configuration file, if specified.") + stop("Error: parameter 'configfile' must be a character string with the path to an s2dverification configuration file, if specified.") } # varmin @@ -1375,6 +1383,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, dims2define <- TRUE is_file_per_member_exp <- rep(nmod, FALSE) exp_work_pieces <- list() + first_time_step_list <- NULL jmod <- 1 while (jmod <= nmod) { first_dataset_file_found <- FALSE @@ -1420,6 +1429,21 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, if (is_file_per_member_exp[jmod]) { replace_values[["MEMBER_NUMBER"]] <- '*' } + if (jsdate == 1) { + work_piecetime <- list(dataset_type = dataset_type, + filename = .ConfigReplaceVariablesInString(quasi_final_path, + replace_values), + namevar = namevar, grid = grid, remap = remap, + remapcells = remapcells, + is_file_per_member = is_file_per_member_exp[jmod], + is_file_per_dataset = FALSE, + lon_limits = c(lonmin, lonmax), + lat_limits = c(latmin, latmax), dimnames = exp[[jmod]][['dimnames']], + single_dataset = single_dataset) + looking_time <- .LoadDataFile(work_piecetime, explore_dims = TRUE, + silent = silent) + first_time_step_list <- c(first_time_step_list, list(looking_time$time_dim)) + } # If the dimensions of the output matrices are still to define, we try to read # the metadata of the data file that corresponds to the current iteration if (dims2define) { @@ -1524,6 +1548,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, jsdate <- jsdate + 1 } replace_values[extra_vars] <- NULL + #first_dataset_file_found <- FALSE jmod <- jmod + 1 } if (dims2define && length(exp) > 0) { @@ -1546,7 +1571,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, # the current date. if (is.null(exp) || dims == 0) { if (is.null(leadtimemax)) { - diff <- Sys.time() - as.POSIXct(sdates[1], format = '%Y%m%d') + diff <- Sys.time() - as.POSIXct(sdates[1], format = '%Y%m%d', tz = "UTC") if (diff > 0) { .warning("Loading observations only and no 'leadtimemax' specified. Data will be loaded from each starting date to current time.") if (storefreq == 'monthly') { @@ -1564,7 +1589,52 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, } leadtimes <- seq(leadtimemin, leadtimemax, sampleperiod) } - + # If there are differences in the first time stamp in exp files: + if (!is.null(exp)) { + in_date <- lapply(first_time_step_list, function(x) { + origin <- as.POSIXct( + paste(strsplit(x$time_units, " ")[[1]][c(3,4)], + collapse = " "), tz = 'UTC') + units <- strsplit(x$time_units, " ")[[1]][1] + if (units == 'hours') { + exp_first_time_step <- as.POSIXct( + x$first_time_step_in_file * + 3600, origin = origin, tz = 'UTC') + } else if (units == 'days') { + exp_first_time_step <- as.POSIXct( + x$first_time_step_in_file * + 86400, origin = origin, tz = 'UTC') + } + day <- as.numeric(format(exp_first_time_step, "%d")) + return(day) + }) + exp_first_time_step <- min(unlist(in_date)) + if (max(unlist(in_date)) > 1) { + leadtimes <- seq(exp_first_time_step, leadtimemax + max(unlist(in_date)) - 1, + sampleperiod) + } + if (leadtimemin > 1 & length(in_date) > 1) { + lags <- lapply(in_date, function(x) {x - in_date[[1]]}) + new_leadtimemin <- lapply(lags, function(x) {leadtimemin - x}) + new_leadtimemax <- lapply(lags, function(x) {leadtimemax - x}) + jmod <- 2 + npieces <- length(exp_work_pieces)/nmod + while (jmod <= nmod) { + jpiece <- 1 + while (jpiece <= npieces) { + exp_work_pieces[[npieces * (jmod - 1) + jpiece]]$leadtimes <- + seq(new_leadtimemin[[jmod]], new_leadtimemax[[jmod]], sampleperiod) + jpiece <- jpiece + 1 + } + jmod <- jmod + 1 + } + } + lag <- 1 - in_date[[1]] + leadtimes <- seq(leadtimemin - lag, leadtimemax #+ max(unlist(in_date)) + lag, + - lag, + sampleperiod) + exp_first_time_step <- leadtimemin - lag + } # Now we start iterating over observations. We try to find the output matrix # dimensions and we build anyway the work pieces corresponding to the observational # data that time-corresponds the experimental data or the time-steps until the @@ -1628,6 +1698,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, found_data <- .LoadDataFile(work_piece, explore_dims = TRUE, silent = silent) found_dims <- found_data$dims var_long_name <- found_data$var_long_name + first_time_step_list <- c(first_time_step_list, list(found_data$time_dim)) units <- found_data$units if (!is.null(found_dims)) { is_2d_var <- found_data$is_2d_var @@ -1702,7 +1773,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, } day <- as.integer(day) startdate <- as.POSIXct(paste(substr(sdate, 1, 4), '-', - substr(sdate, 5, 6), '-', day, ' 12:00:00', sep = '')) + + substr(sdate, 5, 6), '-', day, ' 12:00:00', sep = ''), tz = "UTC") + (leadtimemin - 1) * 86400 year <- as.integer(substr(startdate, 1, 4)) month <- as.integer(substr(startdate, 6, 7)) @@ -1725,7 +1796,18 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, ## This condition must be fulfilled to put all the month time steps ## in the dimension of length nleadtimes. Otherwise it must be cut: #(length(leadtimes) - 1)*sampleperiod + 1 - (jleadtime - 1)*sampleperiod >= days_in_month - day + 1 - obs_file_indices <- seq(day, min(days_in_month, (length(leadtimes) - jleadtime) * sampleperiod + day), sampleperiod) + + ## The first time step in exp could be different from sdate: + if (jleadtime == 1 & !is.null(exp)) { + if (is.null(first_time_step_list[[1]])) { + stop("Check 'time' variable in the experimental files ", + "since not units or first time step have been found.") + } else { + day <- leadtimes[1] + } + } + obs_file_indices <- seq(day, min(days_in_month, + (length(leadtimes) - jleadtime) * sampleperiod + day), sampleperiod) } else { obs_file_indices <- 1 } @@ -1821,7 +1903,8 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, } if (storefreq == 'daily') { - startdate <- startdate + 86400 * sampleperiod * length(obs_file_indices) + startdate <- startdate + 86400 * sampleperiod * + max(obs_file_indices) year <- as.integer(substr(startdate, 1, 4)) month <- as.integer(substr(startdate, 6, 7)) day <- as.integer(substr(startdate, 9, 10)) @@ -2135,49 +2218,54 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, attr(variable, 'daily_agg_cellfun') <- 'none' attr(variable, 'monthly_agg_cellfun') <- 'none' attr(variable, 'verification_time') <- 'none' - - number_ftime <- NULL + + number_ftime <- NULL if (is.null(var_exp)) { mod_data <- NULL - } else { + } else { dim_reorder <- length(dim_exp):1 dim_reorder[2:3] <- dim_reorder[3:2] old_dims <- dim_exp dim_exp <- dim_exp[dim_reorder] - mod_data <- aperm(array(bigmemory::as.matrix(var_exp), dim = old_dims), dim_reorder) + mod_data <- + aperm(array(bigmemory::as.matrix(var_exp), dim = old_dims), dim_reorder) attr(mod_data, 'dimensions') <- names(dim_exp) names(dim(mod_data)) <- names(dim_exp) number_ftime <- dim_exp[["ftime"]] } - + if (is.null(var_obs)) { obs_data <- NULL - } else { + } else { dim_reorder <- length(dim_obs):1 dim_reorder[2:3] <- dim_reorder[3:2] old_dims <- dim_obs dim_obs <- dim_obs[dim_reorder] - obs_data <- aperm(array(bigmemory::as.matrix(var_obs), dim = old_dims), dim_reorder) + obs_data <- + aperm(array(bigmemory::as.matrix(var_obs), dim = old_dims), dim_reorder) attr(obs_data, 'dimensions') <- names(dim_obs) names(dim(obs_data)) <- names(dim_obs) if (is.null(number_ftime)) { number_ftime <- dim_obs[["ftime"]] } } - if (is.null(latitudes)) { lat <- 0 attr(lat, 'cdo_grid_name') <- 'none' attr(lat, 'first_lat') <- 'none' - attr(lat, 'last_lat') <- 'none' + attr(lat, 'last_lat') <- 'none' } else { lat <- latitudes - attr(lat, 'cdo_grid_name') <- if (is.null(grid)) 'none' else grid + attr(lat, 'cdo_grid_name') <- + if (is.null(grid)) + 'none' + else + grid attr(lat, 'first_lat') <- tail(lat, 1) attr(lat, 'last_lat') <- head(lat, 1) } attr(lat, 'projection') <- 'none' - + if (is.null(longitudes)) { lon <- 0 attr(lon, 'cdo_grid_name') <- 'none' @@ -2187,14 +2275,18 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, attr(lon, 'last_lon') <- 'none' } else { lon <- longitudes - attr(lon, 'cdo_grid_name') <- if (is.null(grid)) 'none' else grid + attr(lon, 'cdo_grid_name') <- + if (is.null(grid)) + 'none' + else + grid attr(lon, 'data_across_gw') <- data_across_gw attr(lon, 'array_across_gw') <- array_across_gw attr(lon, 'first_lon') <- lon[which.min(abs(lon - lonmin))] attr(lon, 'last_lon') <- lon[which.min(abs(lon - lonmax))] } attr(lon, 'projection') <- 'none' - + dates <- list() ## we must put a start and end time for each prediction c(start date, forecast time) if (storefreq == 'minutely') { @@ -2206,22 +2298,42 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, } else if (storefreq == 'monthly') { store_period <- 'month' } - + addTime <- function(date, period, n = 1) { seq(date, by = paste(n, period), length = 2)[2] } - + # We build dates, a list with components start and end. # Start is a list with as many components as start dates. # Each component is a vector of the initial POSIXct date of each # forecast time step - dates[["start"]] <- do.call(c, lapply(sdates, + if (!is.null(exp)) { + if (storefreq == 'daily' & leadtimes[[1]] > 1) { + origin <- leadtimes[[1]] - 1 + leadtimemin <- 1 + } else { + origin <- 0 + } + dates[["start"]] <- do.call(c, lapply(sdates, + function(x) { + do.call(c, lapply((origin:(origin + number_ftime - 1)) * sampleperiod, + function(y) { + addTime(as.POSIXct(x, format = "%Y%m%d", tz = "UTC"), + store_period, y + leadtimemin - 1) + })) + })) + } else { + origin <- 0 + dates[["start"]] <- do.call(c, lapply(sdates, function(x) { do.call(c, lapply((0:(number_ftime - 1)) * sampleperiod, function(y) { - addTime(as.POSIXct(x, format = "%Y%m%d"), store_period, y + leadtimemin - 1) + addTime(as.POSIXct(x, format = "%Y%m%d", tz = "UTC"), + store_period, y + leadtimemin - 1) })) })) + } + attr(dates[["start"]], "tzone") <- "UTC" # end is similar to start, but contains the end dates of each forecast # time step dates[["end"]] <- do.call(c, lapply(dates[["start"]], @@ -2231,6 +2343,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, addTime(y, store_period) })) })) + attr(dates[["end"]], "tzone") <- "UTC" tags_to_find <- c('START_DATE', 'MEMBER_NUMBER', 'YEAR', 'MONTH', 'DAY') position_of_tags <- na.omit(match(tags_to_find, names(replace_values))) @@ -2242,14 +2355,19 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, models <- list() for (jmod in 1:length(exp)) { member_names <- paste0("Member_", 1:nmember[jmod]) - models[[exp[[jmod]][["name"]]]] <- list( - InitializationDates = lapply(member_names, + + InitializationDates <- do.call(c, lapply(member_names, function(x) { do.call(c, lapply(sdates, function(y) { - as.POSIXct(y, format = "%Y%m%d") + as.POSIXct(y, format = "%Y%m%d", tz = "UTC") })) - }), + })) + attr(InitializationDates, "tzone") <- "UTC" + + models[[exp[[jmod]][["name"]]]] <- list( + InitializationDates = InitializationDates, Members = member_names) + names(models[[exp[[jmod]][["name"]]]]$InitializationDates) <- member_names attr(models[[exp[[jmod]][["name"]]]], 'dataset') <- exp[[jmod]][["name"]] attr(models[[exp[[jmod]][["name"]]]], 'source') <- { diff --git a/R/Utils.R b/R/Utils.R index 86b2644..1dbf487 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -1,8 +1,3 @@ -#'@importFrom abind abind -#'@importFrom plyr take -#'@importFrom grDevices png jpeg pdf svg bmp tiff -#'@import ncdf4 - ## Function to tell if a regexpr() match is a complete match to a specified name .IsFullMatch <- function(x, name) { ifelse(x > 0 && attributes(x)$match.length == nchar(name), TRUE, FALSE) @@ -211,7 +206,11 @@ if (Sys.which("cdo")[[1]] == "") { stop("Error: CDO libraries not available") } - cdo_version <- as.numeric_version(strsplit(suppressWarnings(system2("cdo", args = '-V', stderr = TRUE))[[1]], ' ')[[1]][5]) + + cdo_version <- strsplit(suppressWarnings(system2("cdo", args = '-V', stderr = TRUE))[[1]], ' ')[[1]][5] + + cdo_version <- as.numeric_version(unlist(strsplit(cdo_version, "[A-Za-z]", fixed = FALSE))[[1]]) + } # If the variable to load is 2-d, we need to determine whether: # - interpolation is needed @@ -596,6 +595,9 @@ nltime <- fnc$var[[namevar]][['dim']][[match(time_dimname, var_dimnames)]]$len expected_dims <- c(expected_dims, time_dimname) dim_matches <- match(expected_dims, var_dimnames) + first_time_step_in_file <- fnc$var[[namevar]][['dim']][[match(time_dimname, + var_dimnames)]]$vals[1] + time_units <- fnc$var[[namevar]][['dim']][[match(time_dimname, var_dimnames)]]$units } else { if (!is.null(old_members_dimname)) { expected_dims[which(expected_dims == 'lev')] <- old_members_dimname @@ -636,7 +638,7 @@ time_indices <- ts(time_indices, start = c(years[1], mons[1]), end = c(years[length(years)], mons[length(mons)]), frequency = 12) - ltimes_list <- list() + ltimes_list <- list() for (sdate in work_piece[['startdates']]) { selected_time_indices <- window(time_indices, start = c(as.numeric( substr(sdate, 1, 4)), as.numeric(substr(sdate, 5, 6))), @@ -881,11 +883,11 @@ weights <- weights / mean(weights, na.rm = TRUE) mean(x * weights, na.rm = TRUE) } else { - weights <- weights / InsertDim(MeanDims(weights, 2, na.rm = TRUE), 2, length(final_lats)) - MeanDims(x * weights, 2, na.rm = TRUE) + weights <- weights / InsertDim(Mean1Dim(weights, 2, narm = TRUE), 2, length(final_lats)) + Mean1Dim(x * weights, 2, narm = TRUE) } } else if (output == 'lat') { - MeanDims(x, 1, na.rm = TRUE) + Mean1Dim(x, 1, narm = TRUE) } else if (output == 'lonlat') { signif(x, 5) } @@ -943,7 +945,9 @@ if (explore_dims) { list(dims = dims, is_2d_var = is_2d_var, grid = grid_name, units = units, var_long_name = var_long_name, - data_across_gw = data_across_gw, array_across_gw = array_across_gw) + data_across_gw = data_across_gw, array_across_gw = array_across_gw, + time_dim = list(first_time_step_in_file = first_time_step_in_file, + time_units = time_units)) } else { ###if (!silent && !is.null(progress_connection) && !is.null(work_piece[['progress_amount']])) { ### foobar <- writeBin(work_piece[['progress_amount']], progress_connection) @@ -977,7 +981,7 @@ lead_times_position <- 4 if (output == 'lonlat') { - sampleData <- s2dv::sampleMap + sampleData <- s2dverification::sampleMap if (is.null(leadtimemax)) { leadtimemax <- dim(sampleData$mod)[lead_times_position] } @@ -988,7 +992,7 @@ dataOut$obs <- sampleData$obs[, , selected_start_dates, selected_lead_times, , ] } else if (output == 'areave') { - sampleData <- s2dv::sampleTimeSeries + sampleData <- s2dverification::sampleTimeSeries if (is.null(leadtimemax)) { leadtimemax <- dim(sampleData$mod)[lead_times_position] } @@ -1629,34 +1633,3 @@ } array1 } - -# only can be used in Trend(). Needs generalization or be replaced by other function. -.reorder <- function(output, time_dim, dim_names) { - # Add dim name back - if (is.null(dim(output))) { - dim(output) <- c(stats = length(output)) - } else { #is an array - if (length(dim(output)) == 1) { - if (!is.null(names(dim(output)))) { - dim(output) <- c(1, dim(output)) - names(dim(output))[1] <- time_dim - } else { - names(dim(output)) <- time_dim - } - } else { # more than one dim - if (names(dim(output))[1] != "") { - dim(output) <- c(1, dim(output)) - names(dim(output))[1] <- time_dim - } else { #regular case - names(dim(output))[1] <- time_dim - } - } - } - # reorder - pos <- match(dim_names, names(dim(output))) - output <- aperm(output, pos) - names(dim(output)) <- dim_names - names(dim(output))[names(dim(output)) == time_dim] <- 'stats' - return(output) -} - -- GitLab From a44c16345fabd8836001fc6c73fdfbe8f5efbb1a Mon Sep 17 00:00:00 2001 From: aho Date: Mon, 2 Nov 2020 12:17:33 +0100 Subject: [PATCH 2/3] Revise some parts back to s2dv --- R/Load.R | 33 +++++++++++++-------------------- R/Utils.R | 46 +++++++++++++++++++++++++++++++++++++++++----- 2 files changed, 54 insertions(+), 25 deletions(-) diff --git a/R/Load.R b/R/Load.R index 9feec5d..0392c74 100644 --- a/R/Load.R +++ b/R/Load.R @@ -42,7 +42,7 @@ #'the individual grid of each dataset but can also be averaged after #'interpolating into a common grid. See parameters 'grid' and 'method'.\cr #'Once the two arrays are filled by calling this function, other functions in -#'the s2dverification package that receive as inputs data formatted in this +#'the s2dv package that receive as inputs data formatted in this #'data structure can be executed (e.g: \code{Clim()} to compute climatologies, #'\code{Ano()} to compute anomalies, ...).\cr\cr #'Load() has many additional parameters to disable values and trim dimensions @@ -455,7 +455,7 @@ #' Warning: list() compulsory even if loading 1 experimental dataset only!\cr #' E.g., list(array(1, dim = c(num_lons, num_lats))) #'@param maskobs See help on parameter 'maskmod'. -#'@param configfile Path to the s2dverification configuration file from which +#'@param configfile Path to the s2dv configuration file from which #' to retrieve information on location in file system (and other) of datasets.\cr #' If not specified, the configuration file used at BSC-ES will be used #' (it is included in the package).\cr @@ -670,13 +670,6 @@ #' } #' } #' -#'@keywords datagen -#'@author History:\cr -#'0.1 - 2011-03 (V. Guemas, \email{virginie.guemas@bsc.es}) - Original code\cr -#'1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens@bsc.es}) - Formatting to CRAN\cr -#'1.2 - 2015-02 (N. Manubens, \email{nicolau.manubens@bsc.es}) - Generalisation + parallelisation\cr -#'1.3 - 2015-07 (N. Manubens, \email{nicolau.manubens@bsc.es}) - Improvements related to configuration file mechanism\cr -#'1.4 - 2016-01 (N. Manubens, \email{nicolau.manubens@bsc.es}) - Added subsetting capabilities #'@examples #'# Let's assume we want to perform verification with data of a variable #'# called 'tos' from a model called 'model' and observed data coming from @@ -762,7 +755,7 @@ #'# Example 1: Providing lists of lists to 'exp' and 'obs': #'# #' \dontrun{ -#'data_path <- system.file('sample_data', package = 's2dverification') +#'data_path <- system.file('sample_data', package = 's2dv') #'exp <- list( #' name = 'experiment', #' path = file.path(data_path, 'model/$EXP_NAME$/monthly_mean', @@ -788,7 +781,7 @@ #'# writing a configuration file). #'# #' \dontrun{ -#'data_path <- system.file('sample_data', package = 's2dverification') +#'data_path <- system.file('sample_data', package = 's2dv') #'expA <- list(name = 'experiment', path = file.path(data_path, #' 'model/$EXP_NAME$/$STORE_FREQ$_mean/$VAR_NAME$_3hourly', #' '$VAR_NAME$_$START_DATE$.nc')) @@ -802,7 +795,7 @@ #' output = 'areave', latmin = 27, latmax = 48, #' lonmin = -12, lonmax = 40) #'# -#'# Example 2: providing character strings in 'exp' and 'obs', and providing +#'# Example 3: providing character strings in 'exp' and 'obs', and providing #'# a configuration file. #'# The configuration file 'sample.conf' that we will create in the example #'# has the proper entries to load these (see ?LoadConfigFile for details on @@ -813,7 +806,7 @@ #'c <- ConfigFileOpen(configfile) #'c <- ConfigEditDefinition(c, 'DEFAULT_VAR_MIN', '-1e19', confirm = FALSE) #'c <- ConfigEditDefinition(c, 'DEFAULT_VAR_MAX', '1e19', confirm = FALSE) -#'data_path <- system.file('sample_data', package = 's2dverification') +#'data_path <- system.file('sample_data', package = 's2dv') #'exp_data_path <- paste0(data_path, '/model/$EXP_NAME$/') #'obs_data_path <- paste0(data_path, '/$OBS_NAME$/') #'c <- ConfigAddEntry(c, 'experiments', dataset_name = 'experiment', @@ -832,11 +825,11 @@ #' } #' \dontshow{ #'startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101') -#'sampleData <- s2dverification:::.LoadSampleData('tos', c('experiment'), -#' c('observation'), startDates, -#' output = 'areave', -#' latmin = 27, latmax = 48, -#' lonmin = -12, lonmax = 40) +#'sampleData <- s2dv:::.LoadSampleData('tos', c('experiment'), +#' c('observation'), startDates, +#' output = 'areave', +#' latmin = 27, latmax = 48, +#' lonmin = -12, lonmax = 40) #' } #'@import parallel bigmemory methods #'@importFrom stats ts window na.omit @@ -1171,9 +1164,9 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, # configfile if (is.null(configfile)) { - configfile <- system.file("config", "BSC.conf", package = "s2dverification") + configfile <- system.file("config", "BSC.conf", package = "s2dv") } else if (!is.character(configfile) || !(nchar(configfile) > 0)) { - stop("Error: parameter 'configfile' must be a character string with the path to an s2dverification configuration file, if specified.") + stop("Error: parameter 'configfile' must be a character string with the path to an s2dv configuration file, if specified.") } # varmin diff --git a/R/Utils.R b/R/Utils.R index 1dbf487..f7dfddb 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -1,3 +1,8 @@ +#'@importFrom abind abind +#'@importFrom plyr take +#'@importFrom grDevices png jpeg pdf svg bmp tiff +#'@import ncdf4 + ## Function to tell if a regexpr() match is a complete match to a specified name .IsFullMatch <- function(x, name) { ifelse(x > 0 && attributes(x)$match.length == nchar(name), TRUE, FALSE) @@ -883,11 +888,11 @@ weights <- weights / mean(weights, na.rm = TRUE) mean(x * weights, na.rm = TRUE) } else { - weights <- weights / InsertDim(Mean1Dim(weights, 2, narm = TRUE), 2, length(final_lats)) - Mean1Dim(x * weights, 2, narm = TRUE) + weights <- weights / InsertDim(MeanDims(weights, 2, narm = TRUE), 2, length(final_lats)) + MeanDims(x * weights, 2, narm = TRUE) } } else if (output == 'lat') { - Mean1Dim(x, 1, narm = TRUE) + MeanDims(x, 1, narm = TRUE) } else if (output == 'lonlat') { signif(x, 5) } @@ -981,7 +986,7 @@ lead_times_position <- 4 if (output == 'lonlat') { - sampleData <- s2dverification::sampleMap + sampleData <- s2dv::sampleMap if (is.null(leadtimemax)) { leadtimemax <- dim(sampleData$mod)[lead_times_position] } @@ -992,7 +997,7 @@ dataOut$obs <- sampleData$obs[, , selected_start_dates, selected_lead_times, , ] } else if (output == 'areave') { - sampleData <- s2dverification::sampleTimeSeries + sampleData <- s2dv::sampleTimeSeries if (is.null(leadtimemax)) { leadtimemax <- dim(sampleData$mod)[lead_times_position] } @@ -1633,3 +1638,34 @@ } array1 } + +# only can be used in Trend(). Needs generalization or be replaced by other function. +.reorder <- function(output, time_dim, dim_names) { + # Add dim name back + if (is.null(dim(output))) { + dim(output) <- c(stats = length(output)) + } else { #is an array + if (length(dim(output)) == 1) { + if (!is.null(names(dim(output)))) { + dim(output) <- c(1, dim(output)) + names(dim(output))[1] <- time_dim + } else { + names(dim(output)) <- time_dim + } + } else { # more than one dim + if (names(dim(output))[1] != "") { + dim(output) <- c(1, dim(output)) + names(dim(output))[1] <- time_dim + } else { #regular case + names(dim(output))[1] <- time_dim + } + } + } + # reorder + pos <- match(dim_names, names(dim(output))) + output <- aperm(output, pos) + names(dim(output)) <- dim_names + names(dim(output))[names(dim(output)) == time_dim] <- 'stats' + return(output) +} + -- GitLab From ad193be584d4882048e3f028380e32471e208542 Mon Sep 17 00:00:00 2001 From: aho Date: Mon, 2 Nov 2020 12:22:30 +0100 Subject: [PATCH 3/3] Typo fix for MeanDims() in Utils --- R/Utils.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/Utils.R b/R/Utils.R index f7dfddb..6e8e713 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -888,11 +888,11 @@ weights <- weights / mean(weights, na.rm = TRUE) mean(x * weights, na.rm = TRUE) } else { - weights <- weights / InsertDim(MeanDims(weights, 2, narm = TRUE), 2, length(final_lats)) - MeanDims(x * weights, 2, narm = TRUE) + weights <- weights / InsertDim(MeanDims(weights, 2, na.rm = TRUE), 2, length(final_lats)) + MeanDims(x * weights, 2, na.rm = TRUE) } } else if (output == 'lat') { - MeanDims(x, 1, narm = TRUE) + MeanDims(x, 1, na.rm = TRUE) } else if (output == 'lonlat') { signif(x, 5) } -- GitLab