From 49556fd381ad364c37106818dc698ecc289cc752 Mon Sep 17 00:00:00 2001 From: Victoria Agudetse Roures Date: Wed, 28 Jun 2023 13:18:30 +0200 Subject: [PATCH 1/3] Change module function names --- modules/Loading/R/load_decadal.R | 554 ++++++++++++++++++++++++++++++ modules/Loading/R/load_seasonal.R | 452 ++++++++++++++++++++++++ 2 files changed, 1006 insertions(+) create mode 100644 modules/Loading/R/load_decadal.R create mode 100644 modules/Loading/R/load_seasonal.R diff --git a/modules/Loading/R/load_decadal.R b/modules/Loading/R/load_decadal.R new file mode 100644 index 00000000..7d0b2d8b --- /dev/null +++ b/modules/Loading/R/load_decadal.R @@ -0,0 +1,554 @@ +# Loading module: +# 1. archive.yml +# 2. recipe.yml +# 3. Load_decadal.R (V) +#setwd('/esarchive/scratch/aho/git/auto-s2s/') + +## TODO: remove paths to personal scratchs +source("/esarchive/scratch/vagudets/repos/csoperational/R/get_regrid_params.R") +# Load required libraries/funs +source("modules/Loading/helper_loading_decadal.R") +source("modules/Loading/R/dates2load.R") +source("modules/Loading/R/check_latlon.R") +source("modules/Loading/R/get_timeidx.R") +source("tools/libs.R") + +#==================================================================== + +# recipe_file <- "recipes/atomic_recipes/recipe_decadal.yml" +# recipe_file <- "recipes/atomic_recipes/recipe_decadal_daily.yml" + +load_decadal <- function(recipe) { + + ## + archive <- read_yaml(paste0("conf/archive_decadal.yml"))$esarchive + + # Print Start() info or not + DEBUG <- FALSE + + ## TODO: this should come from the main script + # Create output folder and log: + + #------------------------- + # Read from recipe: + #------------------------- + exp.name <- recipe$Analysis$Datasets$System$name #'HadGEM3' + ref.name <- recipe$Analysis$Datasets$Reference$name #'era5' + member <- strsplit(recipe$Analysis$Datasets$System$member, ',')[[1]] #c("r1i1p1f2", "r2i1p1f2") +# variable <- recipe$Analysis$Variables$name #'tas' + variable <- strsplit(recipe$Analysis$Variables$name, ", | |,")[[1]] + store.freq <- recipe$Analysis$Variables$freq #monthly_mean + lats.min <- as.numeric(recipe$Analysis$Region$latmin) #0 + lats.max <- as.numeric(recipe$Analysis$Region$latmax) #10 + lons.min <- as.numeric(recipe$Analysis$Region$lonmin) #0 + lons.max <- as.numeric(recipe$Analysis$Region$lonmax) #10 + + # change to: sdates <- dates2load(recipe, logger) + sdates_hcst <- as.numeric(recipe$Analysis$Time$hcst_start):as.numeric(recipe$Analysis$Time$hcst_end) #1960:2015 + sdates_fcst <- recipe$Analysis$Time$fcst + + if (store.freq == "monthly_mean") { + time_ind <- (as.numeric(recipe$Analysis$Time$ftime_min):as.numeric(recipe$Analysis$Time$ftime_max)) + + } else if (store.freq == "daily_mean") { + time_ind <- get_daily_time_ind(ftimemin = as.numeric(recipe$Analysis$Time$ftime_min), + ftimemax = as.numeric(recipe$Analysis$Time$ftime_max), + initial_month = archive$System[[exp.name]]$initial_month, + sdates = sdates_hcst, + calendar = archive$System[[exp.name]]$calendar) + } + +#NOTE: May be used in the future +# season <- recipe$Analysis$Time$season + + #------------------------- + # Read from archive: + #------------------------- + if (store.freq == "monthly_mean") { + table <- archive$System[[exp.name]][[store.freq]]$table[variable] #list(tas = 'Amon') + } else { + table <- 'day' + } + grid <- archive$System[[exp.name]][[store.freq]]$grid[variable] #list(tas = 'gr') + version <- archive$System[[exp.name]][[store.freq]]$version[variable] #list(tas = 'v20210910') + if (identical(member, 'all')) { + member <- strsplit(archive$System[[exp.name]]$member, ',')[[1]] + } + + #------------------------- + # derived from above: + #------------------------- + # Check lat and lon and decide CircularSort + circularsort <- check_latlon(latmin = lats.min, latmax = lats.max, lonmin = lons.min, lonmax = lons.max) + + # generate transform params for system and ref + regrid_params <- get_regrid_params(recipe, archive) + + # Only if the time length in each chunk may differ that we need largest_dims_length to be TRUE. Otherwise, set FALSE to increase efficiency. + need_largest_dims_length <- ifelse(exp.name == 'EC-Earth3-i2', TRUE, FALSE) + + + #------------------------------------------- + # Step 1: Load the hcst + #------------------------------------------- + #monthly and daily + tmp <- get_dcpp_path(archive = archive, exp.name = exp.name, table = table, grid = grid, + version = version, sdates = sdates_hcst) + path_list <- tmp$path_list + multi_path <- tmp$multi_path + + #TODO: to make this case work; enhance Start() if it's possible + if (multi_path & length(variable) > 1) { + stop("The recipe requests multiple variables and start dates from both dpccA-hindcast and dcppB-forecast. This case is not available for now.") + } + + Start_default_arg_list <- list( + dat = path_list, + var = variable, + syear = paste0(sdates_hcst), + chunk = 'all', + chunk_depends = 'syear', + time = indices(time_ind), + time_across = 'chunk', + merge_across_dims = TRUE, + largest_dims_length = need_largest_dims_length, + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(decreasing = TRUE), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = circularsort, + ensemble = member, + transform = regrid_params$fcst.transform, + transform_extra_cells = 2, + transform_params = list(grid = regrid_params$fcst.gridtype, + method = regrid_params$fcst.gridmethod), + transform_vars = c('latitude', 'longitude'), +# path_glob_permissive = 2, # for version + synonims = list(longitude = c('lon', 'longitude'), + latitude = c('lat', 'latitude')), + return_vars = list(latitude = NULL, longitude = NULL, + time = c('syear', 'chunk')), + silent = !DEBUG, + retrieve = T) + + if (length(variable) > 1) { + Start_default_arg_list <- c(Start_default_arg_list, + list(table = table, grid = grid, version = version, + table_depends = 'var', grid_depends = 'var', version_depends = 'var', + metadata_dims = 'var')) + } + + if (!multi_path) { + Start_hcst_arg_list <- Start_default_arg_list + hcst <- do.call(Start, Start_hcst_arg_list) + + } else { + Start_hcst_arg_list <- Start_default_arg_list + Start_hcst_arg_list[['syear']] <- NULL + Start_hcst_arg_list[['chunk_depends']] <- NULL + remove_ind <- which(Start_hcst_arg_list[['return_vars']][['time']] == 'syear') + Start_hcst_arg_list[['return_vars']][['time']] <- Start_hcst_arg_list[['return_vars']][['time']][-remove_ind] + + hcst <- do.call(Start, Start_hcst_arg_list) + + # Reshape and reorder dimensions + ## dat should be 1, syear should be length of dat; reorder dimensions + dim(hcst) <- c(dat = 1, syear = as.numeric(dim(hcst))[1], dim(hcst)[2:6]) + hcst <- s2dv::Reorder(hcst, c('dat', 'var', 'syear', 'time', 'latitude', 'longitude', 'ensemble')) + + # Manipulate time attr because Start() cannot read it correctly + wrong_time_attr <- attr(hcst, 'Variables')$common$time # dim: [time], the first syear only + tmp <- array(dim = c(dim(hcst)[c('syear', 'time')])) + tmp[1, ] <- wrong_time_attr + yr_diff <- (sdates_hcst - sdates_hcst[1])[-1] #diff(sdates_hcst) + for (i_syear in 1:length(yr_diff)) { + tmp[(i_syear + 1), ] <- wrong_time_attr + lubridate::years(yr_diff[i_syear]) + } + attr(hcst, 'Variables')$common$time <- as.POSIXct(tmp, origin = '1970-01-01', tz = 'UTC') + + } + + tmp_time_attr <- attr(hcst, 'Variables')$common$time + + # change syear to c(sday, sweek, syear) + # dim(hcst) should be [dat, var, sday, sweek, syear, time, latitude, longitude, ensemble] + dim(hcst) <- c(dim(hcst)[1:2], sday = 1, sweek = 1, dim(hcst)[3:7]) + if (!identical(dim(tmp_time_attr), dim(hcst)[c('syear', 'time')])) { + error(recipe$Run$logger, + "hcst has problem in matching data and time attr dimension.") + stop() + } + dim(attr(hcst, 'Variables')$common$time) <- c(sday = 1, sweek = 1, dim(tmp_time_attr)) + + #TODO: as.s2dv_cube() needs to be improved to recognize "variable" is under $dat1 + if (multi_path) { + attributes(hcst)$Variables$common[[variable]] <- attributes(hcst)$Variables$dat1[[variable]] + } + + # Change class from startR_array to s2dv_cube + suppressWarnings( + hcst <- as.s2dv_cube(hcst) + ) + +#------------------------------------------- +# Step 2: Load the fcst +#------------------------------------------- + if (!is.null(recipe$Analysis$Time$fcst)) { + + tmp <- get_dcpp_path(archive = archive, exp.name = exp.name, table = table, grid = grid, + version = version, sdates = sdates_fcst) + path_list <- tmp$path_list + multi_path <- tmp$multi_path + + #TODO: to make this case work; enhance Start() if it's possible + if (multi_path & length(variable) > 1) { + stop("The recipe requests multiple variables and start dates from both dpccA-hindcast and dcppB-forecast. This case is not available for now.") + } + + # monthly & daily + if (!multi_path) { + #NOTE: the adjustment for two cases (multiple files per sdate or not) has been made in hcst + Start_fcst_arg_list <- Start_default_arg_list + Start_fcst_arg_list[['dat']] <- path_list + Start_fcst_arg_list[['syear']] <- paste0(sdates_fcst) + fcst <- do.call(Start, Start_fcst_arg_list) + + + } else { # multi_path + + #TODO: time attribute is not correct. Improve Start(). + Start_fcst_arg_list <- Start_default_arg_list + Start_fcst_arg_list[['dat']] <- path_list + Start_fcst_arg_list[['syear']] <- NULL + Start_fcst_arg_list[['chunk_depends']] <- NULL + remove_ind <- which(Start_fcst_arg_list[['return_vars']][['time']] == 'syear') + Start_fcst_arg_list[['return_vars']][['time']] <- Start_fcst_arg_list[['return_vars']][['time']][-remove_ind] + fcst <- do.call(Start, Start_fcst_arg_list) + + # Reshape and reorder dimensions + ## dat should be 1, syear should be length of dat; reorder dimensions + ## dim(fcst) should be [dat, var, syear, time, latitude, longitude, ensemble] + dim(fcst) <- c(dat = 1, syear = as.numeric(dim(fcst))[1], dim(fcst)[2:6]) + fcst <- s2dv::Reorder(fcst, c('dat', 'var', 'syear', 'time', 'latitude', 'longitude', 'ensemble')) + + # Manipulate time attr because Start() cannot read it correctly + wrong_time_attr <- attr(fcst, 'Variables')$common$time # dim: [time], the first syear only + tmp <- array(dim = c(dim(fcst)[c('syear', 'time')])) + tmp[1, ] <- wrong_time_attr + yr_diff <- (sdates_fcst - sdates_fcst[1])[-1] #diff(sdates_fcst) + for (i_syear in 1:length(yr_diff)) { + tmp[(i_syear + 1), ] <- wrong_time_attr + lubridate::years(yr_diff[i_syear]) + } + attr(fcst, 'Variables')$common$time <- as.POSIXct(tmp, origin = '1970-01-01', tz = 'UTC') + + } + + tmp_time_attr <- attr(fcst, 'Variables')$common$time + + # change syear to c(sday, sweek, syear) + # dim(fcst) should be [dat, var, sday, sweek, syear, time, latitude, longitude, ensemble] + dim(fcst) <- c(dim(fcst)[1:2], sday = 1, sweek = 1, dim(fcst)[3:7]) + if (!identical(dim(tmp_time_attr), dim(fcst)[c('syear', 'time')])) { + error(recipe$Run$logger, + "fcst has problem in matching data and time attr dimension.") + stop() + } + dim(attr(fcst, 'Variables')$common$time) <- c(sday = 1, sweek = 1, dim(tmp_time_attr)) + + #TODO: as.s2dv_cube() needs to be improved to recognize "variable" is under $dat1 + if (multi_path) { + attributes(fcst)$Variables$common[[variable]] <- attributes(fcst)$Variables$dat1[[variable]] + } + + # Change class from startR_array to s2dv_cube + suppressWarnings( + fcst <- as.s2dv_cube(fcst) + ) + + # Only syear could be different + if (!identical(dim(hcst$data)[-5], dim(fcst$data)[-5])) { + error(recipe$Run$logger, + "hcst and fcst do not share the same dimension structure.") + stop() + } + + } else { + fcst <- NULL + } + +#------------------------------------------- +# Step 3. Load the reference +#------------------------------------------- + obs.path <- file.path(archive$src, archive$Reference[[ref.name]]$src, + store.freq, "$var$$var_dir$", "$var$_$file_date$.nc") + var_dir_obs <- archive$Reference[[ref.name]][[store.freq]][variable] # list(tas = "_f1h-r1440x721cds", tos = "_f1h-r1440x721cds") + +# obs.path <- file.path(archive$src, archive$Reference[[ref.name]]$src, store.freq, +# paste0(variable, archive$Reference[[ref.name]][[store.freq]][[variable]])) +# obs.files <- paste0('$var$_$file_date$.nc') + + # Get from startR_cube +# dates <- attr(hcst, 'Variables')$common$time + # Get from s2dv_cube + dates <- hcst$attrs$Dates + dates_file <- sapply(dates, format, '%Y%m') + dim(dates_file) <- dim(dates) + + if (store.freq == "daily_mean") { +#//////////////// +# Method 1: use hcst time attr as obs time selector +#//////////////// + + # Set hour to 12:00 to ensure correct date retrieval for daily data + lubridate::hour(dates) <- 12 + lubridate::minute(dates) <- 00 + # Restore correct dimensions + dim(dates) <- dim(dates_file) + + obs <- Start(dat = obs.path, + var = variable, + var_dir = var_dir_obs, + var_dir_depends = 'var', + file_date = unique(format(dates, '%Y%m')), + time = dates, # [sday, sweek, syear, time] + time_across = 'file_date', + merge_across_dims = TRUE, + split_multiselected_dims = TRUE, + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(decreasing = TRUE), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = circularsort, + transform = regrid_params$obs.transform, + transform_extra_cells = 2, + transform_params = list(grid = regrid_params$obs.gridtype, #nc file + method = regrid_params$obs.gridmethod), + transform_vars = c('latitude', 'longitude'), + synonims = list(latitude = c('lat','latitude'), + longitude = c('lon','longitude')), + return_vars = list(latitude = NULL, longitude = NULL, + time = 'file_date'), + silent = !DEBUG, + retrieve = TRUE) + + } else if (store.freq == "monthly_mean") { +#//////////////// +# Method 2: reshape hcst time attr's date into an array with time dim then as obs date selector +#//////////////// + + obs <- Start(dat = obs.path, + var = variable, + var_dir = var_dir_obs, + var_dir_depends = 'var', + file_date = dates_file, #dates_arr, # [sday, sweek, syear, time] + split_multiselected_dims = TRUE, + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(decreasing = TRUE), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = circularsort, + transform = regrid_params$obs.transform, + transform_extra_cells = 2, + transform_params = list(grid = regrid_params$obs.gridtype, #nc file + method = regrid_params$obs.gridmethod), + transform_vars = c('latitude', 'longitude'), + synonims = list(latitude = c('lat','latitude'), + longitude = c('lon','longitude')), + return_vars = list(latitude = NULL, longitude = NULL, + time = 'file_date'), + metadata_dims = 'var', + silent = !DEBUG, + retrieve = TRUE) + } + + +#dim(attr(obs, 'Variables')$common$time) +# sday sweek syear time +# 1 1 2 14 + + # Remove var_dir dimension + obs <- Subset(obs, along = "var_dir", indices = 1, drop = "selected") + + # Only ensemble dim could be different + if (!identical(dim(obs), dim(hcst$data)[-9])) { + error(recipe$Run$logger, + "obs and hcst dimensions do not match.") + stop() + } + # Add ensemble dim to obs + dim(obs) <- c(dim(obs), ensemble = 1) + + # Change class from startR_array to s2dv_cube + suppressWarnings( + obs <- as.s2dv_cube(obs) + ) + +#------------------------------------------- +# Step 4. Verify the consistance between data +#------------------------------------------- + # dimension + if (any(!names(dim(obs$data)) %in% names(dim(hcst$data)))) { + error(recipe$Run$logger, + "hcst and obs don't share the same dimension names.") + stop() + } else { + ens_ind <- which(names(dim(obs$data)) == 'ensemble') + match_ind <- match(names(dim(obs$data))[-ens_ind], names(dim(hcst$data))) + if (!all(dim(hcst$data)[match_ind] == dim(obs$data)[-ens_ind])) { + error(recipe$Run$logger, + "hcst and obs don't share the same dimension length.") + stop() + } + } + + # time attribute + if (!identical(format(hcst$attrs$Dates, '%Y%m'), + format(obs$attrs$Dates, '%Y%m'))) { + error(recipe$Run$logger, + "hcst and obs don't share the same time.") + stop() + } + + # lat and lon attributes + if (!(recipe$Analysis$Regrid$type == 'none')) { + if (!isTRUE(all.equal(as.vector(hcst$lat), as.vector(obs$lat)))) { + lat_error_msg <- paste("Latitude mismatch between hcst and obs.", + "Please check the original grids and the", + "regrid parameters in your recipe.") + error(recipe$Run$logger, lat_error_msg) + hcst_lat_msg <- paste0("First hcst lat: ", hcst$lat[1], + "; Last hcst lat: ", hcst$lat[length(hcst$lat)]) + info(recipe$Run$logger, hcst_lat_msg) + obs_lat_msg <- paste0("First obs lat: ", obs$lat[1], + "; Last obs lat: ", obs$lat[length(obs$lat)]) + info(recipe$Run$logger, obs_lat_msg) + stop("hcst and obs don't share the same latitudes.") + } + + if (!isTRUE(all.equal(as.vector(hcst$lon), as.vector(obs$lon)))) { + lon_error_msg <- paste("Longitude mismatch between hcst and obs.", + "Please check the original grids and the", + "regrid parameters in your recipe.") + error(recipe$Run$logger, lon_error_msg) + hcst_lon_msg <- paste0("First hcst lon: ", hcst$lon[1], + "; Last hcst lon: ", hcst$lon[length(hcst$lon)]) + info(recipe$Run$logger, hcst_lon_msg) + obs_lon_msg <- paste0("First obs lon: ", obs$lon[1], + "; Last obs lon: ", obs$lon[length(obs$lon)]) + info(recipe$Run$logger, obs_lon_msg) + stop("hcst and obs don't share the same longitudes.") + } + } + + # Check fcst + if (!is.null(fcst)) { + # dimension + if (any(!names(dim(fcst$data)) %in% names(dim(hcst$data)))) { + error(recipe$Run$logger, + "hcst and fcst don't share the same dimension names.") + stop() + } else { + ens_ind <- which(names(dim(fcst$data)) %in% c('ensemble', 'syear')) + match_ind <- match(names(dim(fcst$data))[-ens_ind], names(dim(hcst$data))) + if (!all(dim(hcst$data)[match_ind] == dim(fcst$data)[-ens_ind])) { + error(recipe$Run$logger, + "hcst and fcst don't share the same dimension length.") + stop() + } + } + + # lat and lon attributes + if (!(recipe$Analysis$Regrid$type == 'none')) { + if (!identical(as.vector(hcst$lat), as.vector(fcst$lat))) { + lat_error_msg <- paste("Latitude mismatch between hcst and fcst.", + "Please check the original grids and the", + "regrid parameters in your recipe.") + error(recipe$Run$logger, lat_error_msg) + hcst_lat_msg <- paste0("First hcst lat: ", hcst$lat[1], + "; Last hcst lat: ", hcst$lat[length(hcst$lat)]) + info(recipe$Run$logger, hcst_lat_msg) + fcst_lat_msg <- paste0("First fcst lat: ", fcst$lat[1], + "; Last fcst lat: ", fcst$lat[length(fcst$lat)]) + info(recipe$Run$logger, fcst_lat_msg) + stop("hcst and fcst don't share the same latitudes.") + } + + if (!identical(as.vector(hcst$lon), as.vector(fcst$lon))) { + lon_error_msg <- paste("Longitude mismatch between hcst and fcst.", + "Please check the original grids and the", + "regrid parameters in your recipe.") + error(recipe$Run$logger, lon_error_msg) + hcst_lon_msg <- paste0("First hcst lon: ", hcst$lon[1], + "; Last hcst lon: ", hcst$lon[length(hcst$lon)]) + info(recipe$Run$logger, hcst_lon_msg) + fcst_lon_msg <- paste0("First fcst lon: ", fcst$lon[1], + "; Last fcst lon: ", fcst$lon[length(fcst$lon)]) + info(recipe$Run$logger, fcst_lon_msg) + stop("hcst and fcst don't share the same longitudes.") + } + } + + } + + +#------------------------------------------- +# Step 5. Tune data +#------------------------------------------- + # Remove negative values in accumulative variables + dictionary <- read_yaml("conf/variable-dictionary.yml") + for (var_idx in 1:length(variable)) { + var_name <- variable[var_idx] + if (dictionary$vars[[var_name]]$accum) { + info(recipe$Run$logger, + paste0("Accumulated variable ", var_name, + ": setting negative values to zero.")) + # obs$data[, var_idx, , , , , , , ] <- pmax(Subset(obs$data, + # along = "var", + # indices = var_idx, F), 0) + obs$data[, var_idx, , , , , , , ][obs$data[, var_idx, , , , , , , ] < 0] <- 0 + hcst$data[, var_idx, , , , , , , ][hcst$data[, var_idx, , , , , , , ] < 0] <- 0 + if (!is.null(fcst)) { + fcst$data[, var_idx, , , , , , , ][fcst$data[, var_idx, , , , , , , ] < 0] <- 0 + } + } + + # Convert prlr from m/s to mm/day + ## TODO: Make a unit conversion function + if (variable[[var_idx]] == "prlr") { + # Verify that the units are m/s and the same in obs and hcst + if (((obs$attrs$Variable$metadata[[var_name]]$units == "m s-1") || + (obs$attrs$Variable$metadata[[var_name]]$units == "m s**-1")) && + ((hcst$attrs$Variable$metadata[[var_name]]$units == "m s-1") || + (hcst$attrs$Variable$metadata[[var_name]]$units == "m s**-1"))) { + info(recipe$Run$logger, "Converting precipitation from m/s to mm/day.") + obs$data[, var_idx, , , , , , , ] <- + obs$data[, var_idx, , , , , , , ]*86400*1000 + obs$attrs$Variable$metadata[[var_name]]$units <- "mm/day" + hcst$data[, var_idx, , , , , , , ] <- + hcst$data[, var_idx, , , , , , , ]*86400*1000 + hcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day" + if (!is.null(fcst)) { + fcst$data[, var_idx, , , , , , , ] <- + fcst$data[, var_idx, , , , , , , ]*86400*1000 + fcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day" + } + } + } + } + +#------------------------------------------- +# Step 6. Print summary +#------------------------------------------- + + # Print a summary of the loaded data for the user, for each object + if (recipe$Run$logger$threshold <= 2) { + data_summary(hcst, recipe) + data_summary(obs, recipe) + if (!is.null(fcst)) { + data_summary(fcst, recipe) + } + } + + info(recipe$Run$logger, + "##### DATA LOADING COMPLETED SUCCESSFULLY #####") + + + return(list(hcst = hcst, fcst = fcst, obs = obs)) +} diff --git a/modules/Loading/R/load_seasonal.R b/modules/Loading/R/load_seasonal.R new file mode 100644 index 00000000..c54ef908 --- /dev/null +++ b/modules/Loading/R/load_seasonal.R @@ -0,0 +1,452 @@ +## TODO: remove paths to personal scratchs +source("/esarchive/scratch/vagudets/repos/csoperational/R/get_regrid_params.R") +# Load required libraries/funs +source("modules/Loading/R/dates2load.R") +source("modules/Loading/R/get_timeidx.R") +source("modules/Loading/R/check_latlon.R") +## TODO: Move to prepare_outputs.R +source("tools/libs.R") + +load_seasonal <- function(recipe) { + + # ------------------------------------------- + # Set params ----------------------------------------- + + hcst.inityear <- recipe$Analysis$Time$hcst_start + hcst.endyear <- recipe$Analysis$Time$hcst_end + lats.min <- recipe$Analysis$Region$latmin + lats.max <- recipe$Analysis$Region$latmax + lons.min <- recipe$Analysis$Region$lonmin + lons.max <- recipe$Analysis$Region$lonmax + ref.name <- recipe$Analysis$Datasets$Reference$name + exp.name <- recipe$Analysis$Datasets$System$name + + variable <- strsplit(recipe$Analysis$Variables$name, ", | |,")[[1]] + store.freq <- recipe$Analysis$Variables$freq + + # get sdates array + ## LOGGER: Change dates2load to extract logger from recipe? + sdates <- dates2load(recipe, recipe$Run$logger) + + idxs <- NULL + idxs$hcst <- get_timeidx(sdates$hcst, + recipe$Analysis$Time$ftime_min, + recipe$Analysis$Time$ftime_max, + time_freq=store.freq) + + if (!(is.null(sdates$fcst))) { + idxs$fcst <- get_timeidx(sdates$fcst, + recipe$Analysis$Time$ftime_min, + recipe$Analysis$Time$ftime_max, + time_freq=store.freq) + } + + ## TODO: Examine this verifications part, verify if it's necessary + # stream <- verifications$stream + # sdates <- verifications$fcst.sdate + + ## TODO: define fcst.name + ##fcst.name <- recipe$Analysis$Datasets$System[[sys]]$name + + # get esarchive datasets dict: + ## TODO: Adapt to 'filesystem' option in recipe + archive <- read_yaml("conf/archive.yml")$esarchive + exp_descrip <- archive$System[[exp.name]] + + freq.hcst <- unlist(exp_descrip[[store.freq]][variable[1]]) + reference_descrip <- archive$Reference[[ref.name]] + freq.obs <- unlist(reference_descrip[[store.freq]][variable[1]]) + obs.dir <- reference_descrip$src + fcst.dir <- exp_descrip$src + hcst.dir <- exp_descrip$src + fcst.nmember <- exp_descrip$nmember$fcst + hcst.nmember <- exp_descrip$nmember$hcst + + ## TODO: it is necessary? + ##if ("accum" %in% names(reference_descrip)) { + ## accum <- unlist(reference_descrip$accum[store.freq][[1]]) + ##} else { + ## accum <- FALSE + ##} + + var_dir_obs <- reference_descrip[[store.freq]][variable] + var_dir_exp <- exp_descrip[[store.freq]][variable] + + # ----------- + obs.path <- paste0(archive$src, + obs.dir, store.freq, "/$var$", "$var_dir$", + "/$var$_$file_date$.nc") + + hcst.path <- paste0(archive$src, + hcst.dir, store.freq, "/$var$", "$var_dir$", + "$var$_$file_date$.nc") + + fcst.path <- paste0(archive$src, + hcst.dir, store.freq, "/$var$", "$var_dir$", + "/$var$_$file_date$.nc") + + # Define regrid parameters: + #------------------------------------------------------------------- + regrid_params <- get_regrid_params(recipe, archive) + + # Longitude circular sort and latitude check + #------------------------------------------------------------------- + circularsort <- check_latlon(lats.min, lats.max, lons.min, lons.max) + + if (recipe$Analysis$Variables$freq == "monthly_mean"){ + split_multiselected_dims = TRUE + } else { + split_multiselected_dims = FALSE + } + + # Load hindcast + #------------------------------------------------------------------- + hcst <- Start(dat = hcst.path, + var = variable, + var_dir = var_dir_exp, + file_date = sdates$hcst, + time = idxs$hcst, + var_dir_depends = 'var', + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = circularsort, + transform = regrid_params$fcst.transform, + transform_params = list(grid = regrid_params$fcst.gridtype, + method = regrid_params$fcst.gridmethod), + transform_vars = c('latitude', 'longitude'), + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('lon', 'longitude'), + ensemble = c('member', 'ensemble')), + ensemble = indices(1:hcst.nmember), + metadata_dims = 'var', # change to just 'var'? + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'file_date'), + split_multiselected_dims = split_multiselected_dims, + retrieve = TRUE) + + # Remove var_dir dimension + if ("var_dir" %in% names(dim(hcst))) { + hcst <- Subset(hcst, along = "var_dir", indices = 1, drop = "selected") + } + + if (store.freq %in% c("daily_mean", "daily")) { + # Adjusts dims for daily case, could be removed if startR allows + # multidim split + names(dim(hcst))[which(names(dim(hcst)) == 'file_date')] <- "syear" + default_dims <- c(dat = 1, var = 1, sday = 1, + sweek = 1, syear = 1, time = 1, + latitude = 1, longitude = 1, ensemble = 1) + default_dims[names(dim(hcst))] <- dim(hcst) + dim(hcst) <- default_dims + # Change time attribute dimensions + default_time_dims <- c(sday = 1, sweek = 1, syear = 1, time = 1) + names(dim(attr(hcst, "Variables")$common$time))[which(names( + dim(attr(hcst, "Variables")$common$time)) == 'file_date')] <- "syear" + default_time_dims[names(dim(attr(hcst, "Variables")$common$time))] <- + dim(attr(hcst, "Variables")$common$time) + dim(attr(hcst, "Variables")$common$time) <- default_time_dims + } + + # Convert hcst to s2dv_cube object + ## TODO: Give correct dimensions to $Dates + ## (sday, sweek, syear instead of file_date) + hcst <- as.s2dv_cube(hcst) + # Adjust dates for models where the time stamp goes into the next month + if (recipe$Analysis$Variables$freq == "monthly_mean") { + hcst$attrs$Dates[] <- hcst$attrs$Dates - seconds(exp_descrip$time_stamp_lag) + } + + # Load forecast + #------------------------------------------------------------------- + if (!is.null(recipe$Analysis$Time$fcst_year)) { + # the call uses file_date instead of fcst_syear so that it can work + # with the daily case and the current version of startR not allowing + # multiple dims split + + fcst <- Start(dat = fcst.path, + var = variable, + var_dir = var_dir_exp, + var_dir_depends = 'var', + file_date = sdates$fcst, + time = idxs$fcst, + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = circularsort, + transform = regrid_params$fcst.transform, + transform_params = list(grid = regrid_params$fcst.gridtype, + method = regrid_params$fcst.gridmethod), + transform_vars = c('latitude', 'longitude'), + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('lon', 'longitude'), + ensemble = c('member', 'ensemble')), + ensemble = indices(1:fcst.nmember), + metadata_dims = 'var', + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'file_date'), + split_multiselected_dims = split_multiselected_dims, + retrieve = TRUE) + + if ("var_dir" %in% names(dim(fcst))) { + fcst <- Subset(fcst, along = "var_dir", indices = 1, drop = "selected") + } + + if (store.freq %in% c("daily_mean", "daily")) { + # Adjusts dims for daily case, could be removed if startR allows + # multidim split + names(dim(fcst))[which(names(dim(fcst)) == 'file_date')] <- "syear" + default_dims <- c(dat = 1, var = 1, sday = 1, + sweek = 1, syear = 1, time = 1, + latitude = 1, longitude = 1, ensemble = 1) + default_dims[names(dim(fcst))] <- dim(fcst) + dim(fcst) <- default_dims + # Change time attribute dimensions + default_time_dims <- c(sday = 1, sweek = 1, syear = 1, time = 1) + names(dim(attr(fcst, "Variables")$common$time))[which(names( + dim(attr(fcst, "Variables")$common$time)) == 'file_date')] <- "syear" + default_time_dims[names(dim(attr(fcst, "Variables")$common$time))] <- + dim(attr(fcst, "Variables")$common$time) + dim(attr(fcst, "Variables")$common$time) <- default_time_dims + } + + # Convert fcst to s2dv_cube + fcst <- as.s2dv_cube(fcst) + # Adjust dates for models where the time stamp goes into the next month + if (recipe$Analysis$Variables$freq == "monthly_mean") { + fcst$attrs$Dates[] <- + fcst$attrs$Dates - seconds(exp_descrip$time_stamp_lag) + } + + } else { + fcst <- NULL + } + + # Load reference + #------------------------------------------------------------------- + + # Obtain dates and date dimensions from the loaded hcst data to make sure + # the corresponding observations are loaded correctly. + dates <- hcst$attrs$Dates + dim(dates) <- hcst$dims[c("sday", "sweek", "syear", "time")] + + # Separate Start() call for monthly vs daily data + if (store.freq == "monthly_mean") { + + dates_file <- format(as.Date(dates, '%Y%m%d'), "%Y%m") + dim(dates_file) <- dim(dates) + + obs <- Start(dat = obs.path, + var = variable, + var_dir = var_dir_obs, + var_dir_depends = 'var', + file_date = dates_file, + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = circularsort, + transform = regrid_params$obs.transform, + transform_params = list(grid = regrid_params$obs.gridtype, + method = regrid_params$obs.gridmethod), + transform_vars = c('latitude', 'longitude'), + synonims = list(latitude = c('lat','latitude'), + longitude = c('lon','longitude')), + metadata_dims = 'var', + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'file_date'), + split_multiselected_dims = TRUE, + retrieve = TRUE) + + } else if (store.freq %in% c("daily_mean", "daily")) { + + # Get year and month for file_date + dates_file <- sapply(dates, format, '%Y%m') + dim(dates_file) <- dim(dates) + # Set hour to 12:00 to ensure correct date retrieval for daily data + lubridate::hour(dates) <- 12 + lubridate::minute(dates) <- 00 + # Restore correct dimensions + dim(dates) <- dim(dates_file) + + obs <- Start(dat = obs.path, + var = variable, + var_dir = var_dir_obs, + var_dir_depends = 'var', + file_date = sort(unique(dates_file)), + time = dates, + time_var = 'time', + time_across = 'file_date', + merge_across_dims = TRUE, + merge_across_dims_narm = TRUE, + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = circularsort, + transform = regrid_params$obs.transform, + transform_params = list(grid = regrid_params$obs.gridtype, + method = regrid_params$obs.gridmethod), + transform_vars = c('latitude', 'longitude'), + synonims = list(latitude = c('lat','latitude'), + longitude = c('lon','longitude')), + metadata_dims = 'var', + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'file_date'), + split_multiselected_dims = TRUE, + retrieve = TRUE) + } + + + # Remove var_dir dimension + if ("var_dir" %in% names(dim(obs))) { + obs <- Subset(obs, along = "var_dir", indices = 1, drop = "selected") + } + # Adds ensemble dim to obs (for consistency with hcst/fcst) + default_dims <- c(dat = 1, var = 1, sday = 1, + sweek = 1, syear = 1, time = 1, + latitude = 1, longitude = 1, ensemble = 1) + default_dims[names(dim(obs))] <- dim(obs) + dim(obs) <- default_dims + + # Convert obs to s2dv_cube + obs <- as.s2dv_cube(obs) + + # Check for consistency between hcst and obs grid + if (!(recipe$Analysis$Regrid$type == 'none')) { + if (!isTRUE(all.equal(as.vector(hcst$lat), as.vector(obs$lat)))) { + lat_error_msg <- paste("Latitude mismatch between hcst and obs.", + "Please check the original grids and the", + "regrid parameters in your recipe.") + error(recipe$Run$logger, lat_error_msg) + hcst_lat_msg <- paste0("First hcst lat: ", hcst$lat[1], + "; Last hcst lat: ", hcst$lat[length(hcst$lat)]) + info(recipe$Run$logger, hcst_lat_msg) + obs_lat_msg <- paste0("First obs lat: ", obs$lat[1], + "; Last obs lat: ", obs$lat[length(obs$lat)]) + info(recipe$Run$logger, obs_lat_msg) + stop("hcst and obs don't share the same latitudes.") + } + if (!isTRUE(all.equal(as.vector(hcst$lon), as.vector(obs$lon)))) { + lon_error_msg <- paste("Longitude mismatch between hcst and obs.", + "Please check the original grids and the", + "regrid parameters in your recipe.") + error(recipe$Run$logger, lon_error_msg) + hcst_lon_msg <- paste0("First hcst lon: ", hcst$lon[1], + "; Last hcst lon: ", hcst$lon[length(hcst$lon)]) + info(recipe$Run$logger, hcst_lon_msg) + obs_lon_msg <- paste0("First obs lon: ", obs$lon[1], + "; Last obs lon: ", obs$lon[length(obs$lon)]) + info(recipe$Run$logger, obs_lon_msg) + stop("hcst and obs don't share the same longitudes.") + + } + } + + # Remove negative values in accumulative variables + dictionary <- read_yaml("conf/variable-dictionary.yml") + for (var_idx in 1:length(variable)) { + var_name <- variable[var_idx] + if (dictionary$vars[[var_name]]$accum) { + info(recipe$Run$logger, + paste0("Accumulated variable ", var_name, + ": setting negative values to zero.")) + # obs$data[, var_idx, , , , , , , ] <- pmax(Subset(obs$data, + # along = "var", + # indices = var_idx, F), 0) + obs$data[, var_idx, , , , , , , ][obs$data[, var_idx, , , , , , , ] < 0] <- 0 + hcst$data[, var_idx, , , , , , , ][hcst$data[, var_idx, , , , , , , ] < 0] <- 0 + if (!is.null(fcst)) { + fcst$data[, var_idx, , , , , , , ][fcst$data[, var_idx, , , , , , , ] < 0] <- 0 + } + } + + # Convert prlr from m/s to mm/day + ## TODO: Make a unit conversion function + if (variable[[var_idx]] == "prlr") { + # Verify that the units are m/s and the same in obs and hcst + if (((obs$attrs$Variable$metadata[[var_name]]$units == "m s-1") || + (obs$attrs$Variable$metadata[[var_name]]$units == "m s**-1")) && + ((hcst$attrs$Variable$metadata[[var_name]]$units == "m s-1") || + (hcst$attrs$Variable$metadata[[var_name]]$units == "m s**-1"))) { + info(recipe$Run$logger, "Converting precipitation from m/s to mm/day.") + obs$data[, var_idx, , , , , , , ] <- + obs$data[, var_idx, , , , , , , ]*86400*1000 + obs$attrs$Variable$metadata[[var_name]]$units <- "mm/day" + hcst$data[, var_idx, , , , , , , ] <- + hcst$data[, var_idx, , , , , , , ]*86400*1000 + hcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day" + if (!is.null(fcst)) { + fcst$data[, var_idx, , , , , , , ] <- + fcst$data[, var_idx, , , , , , , ]*86400*1000 + fcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day" + } + } + } + } + # Compute anomalies if requested + # Print a summary of the loaded data for the user, for each object + if (recipe$Run$logger$threshold <= 2) { + data_summary(hcst, recipe) + data_summary(obs, recipe) + if (!is.null(fcst)) { + data_summary(fcst, recipe) + } + } + + info(recipe$Run$logger, + "##### DATA LOADING COMPLETED SUCCESSFULLY #####") + + ############################################################################ + # + # CHECKS ON MISSING FILES + # + ############################################################################ + + #obs.NA_dates.ind <- Apply(obs, + # fun=(function(x){ all(is.na(x))}), + # target_dims=c('time', 'latitude', 'longitude'))[[1]] + #obs.NA_dates <- dates_file[obs.NA_dates.ind] + #obs.NA_dates <- obs.NA_dates[order(obs.NA_dates)] + #obs.NA_files <- paste0(obs.dir, store.freq,"/",variable,"_", + # freq.obs,"obs.grid","/",variable,"_",obs.NA_dates,".nc") + # + #if (any(is.na(hcst))){ + # fatal(recipe$Run$logger, + # paste(" ERROR: MISSING HCST VALUES FOUND DURING LOADING # ", + # " ################################################# ", + # " ###### MISSING FILES #### ", + # " ################################################# ", + # "hcst files:", + # hcst.NA_files, + # " ################################################# ", + # " ################################################# ", + # sep="\n")) + # quit(status = 1) + #} + # + #if (any(is.na(obs)) && !identical(obs.NA_dates,character(0))){ + # fatal(recipe$logger, + # paste(" ERROR: MISSING OBS VALUES FOUND DURING LOADING # ", + # " ################################################# ", + # " ###### MISSING FILES #### ", + # " ################################################# ", + # "obs files:", + # obs.NA_files, + # " ################################################# ", + # " ################################################# ", + # sep="\n")) + # quit(status=1) + #} + # + #info(recipe$logger, + # "######### DATA LOADING COMPLETED SUCCESFULLY ##############") + + ############################################################################ + ############################################################################ + + return(list(hcst = hcst, fcst = fcst, obs = obs)) + +} -- GitLab From cdd16215a0ead1e713a9e704c96b0fbe9a92b5ca Mon Sep 17 00:00:00 2001 From: Victoria Agudetse Roures Date: Wed, 28 Jun 2023 13:19:15 +0200 Subject: [PATCH 2/3] Change module function names (WIP) --- example_scripts/test_seasonal.R | 12 +- modules/Anomalies/Anomalies.R | 2 +- modules/Calibration/Calibration.R | 2 +- modules/Loading/Loading.R | 461 +------------------------- modules/Saving/Saving.R | 11 +- modules/Skill/Skill.R | 4 +- modules/Visualization/Visualization.R | 10 +- 7 files changed, 33 insertions(+), 469 deletions(-) diff --git a/example_scripts/test_seasonal.R b/example_scripts/test_seasonal.R index 1afb9549..22524452 100644 --- a/example_scripts/test_seasonal.R +++ b/example_scripts/test_seasonal.R @@ -9,17 +9,17 @@ recipe_file <- "recipes/atomic_recipes/recipe_test_multivar.yml" recipe <- prepare_outputs(recipe_file) # Load datasets -data <- load_datasets(recipe) +data <- Loading(recipe) # Calibrate datasets -data <- calibrate_datasets(recipe, data) +data <- Calibration(recipe, data) # Compute anomalies -data <- compute_anomalies(recipe, data) +data <- Anomalies(recipe, data) # Compute skill metrics -skill_metrics <- compute_skill_metrics(recipe, data) +skill_metrics <- Skill(recipe, data) # Compute percentiles and probability bins -probabilities <- compute_probabilities(recipe, data) +probabilities <- Probabilities(recipe, data) # Export all data to netCDF ## TODO: Fix plotting # save_data(recipe, data, skill_metrics, probabilities) # Plot data -plot_data(recipe, data, skill_metrics, probabilities, significance = T) +Visualization(recipe, data, skill_metrics, probabilities, significance = T) diff --git a/modules/Anomalies/Anomalies.R b/modules/Anomalies/Anomalies.R index 28092625..969dfebb 100644 --- a/modules/Anomalies/Anomalies.R +++ b/modules/Anomalies/Anomalies.R @@ -1,7 +1,7 @@ # Compute the hcst, obs and fcst anomalies with or without cross-validation # and return them, along with the hcst and obs climatologies. -compute_anomalies <- function(recipe, data) { +Anomalies <- function(recipe, data) { if (is.null(recipe$Analysis$Workflow$Anomalies$compute)) { error(recipe$Run$logger, diff --git a/modules/Calibration/Calibration.R b/modules/Calibration/Calibration.R index dedfbd85..2e0cf7f0 100644 --- a/modules/Calibration/Calibration.R +++ b/modules/Calibration/Calibration.R @@ -1,5 +1,5 @@ -calibrate_datasets <- function(recipe, data) { +Calibration <- function(recipe, data) { # Function that calibrates the hindcast using the method stated in the # recipe. If the forecast is not null, it calibrates it as well. # diff --git a/modules/Loading/Loading.R b/modules/Loading/Loading.R index d945d5af..9b9ee86b 100644 --- a/modules/Loading/Loading.R +++ b/modules/Loading/Loading.R @@ -1,452 +1,15 @@ -## TODO: remove paths to personal scratchs -source("/esarchive/scratch/vagudets/repos/csoperational/R/get_regrid_params.R") -# Load required libraries/funs -source("modules/Loading/R/dates2load.R") -source("modules/Loading/R/get_timeidx.R") -source("modules/Loading/R/check_latlon.R") -## TODO: Move to prepare_outputs.R -source("tools/libs.R") - -load_datasets <- function(recipe) { - - # ------------------------------------------- - # Set params ----------------------------------------- - - hcst.inityear <- recipe$Analysis$Time$hcst_start - hcst.endyear <- recipe$Analysis$Time$hcst_end - lats.min <- recipe$Analysis$Region$latmin - lats.max <- recipe$Analysis$Region$latmax - lons.min <- recipe$Analysis$Region$lonmin - lons.max <- recipe$Analysis$Region$lonmax - ref.name <- recipe$Analysis$Datasets$Reference$name - exp.name <- recipe$Analysis$Datasets$System$name - - variable <- strsplit(recipe$Analysis$Variables$name, ", | |,")[[1]] - store.freq <- recipe$Analysis$Variables$freq - - # get sdates array - ## LOGGER: Change dates2load to extract logger from recipe? - sdates <- dates2load(recipe, recipe$Run$logger) - - idxs <- NULL - idxs$hcst <- get_timeidx(sdates$hcst, - recipe$Analysis$Time$ftime_min, - recipe$Analysis$Time$ftime_max, - time_freq=store.freq) - - if (!(is.null(sdates$fcst))) { - idxs$fcst <- get_timeidx(sdates$fcst, - recipe$Analysis$Time$ftime_min, - recipe$Analysis$Time$ftime_max, - time_freq=store.freq) - } - - ## TODO: Examine this verifications part, verify if it's necessary - # stream <- verifications$stream - # sdates <- verifications$fcst.sdate - - ## TODO: define fcst.name - ##fcst.name <- recipe$Analysis$Datasets$System[[sys]]$name - - # get esarchive datasets dict: - ## TODO: Adapt to 'filesystem' option in recipe - archive <- read_yaml("conf/archive.yml")$esarchive - exp_descrip <- archive$System[[exp.name]] - - freq.hcst <- unlist(exp_descrip[[store.freq]][variable[1]]) - reference_descrip <- archive$Reference[[ref.name]] - freq.obs <- unlist(reference_descrip[[store.freq]][variable[1]]) - obs.dir <- reference_descrip$src - fcst.dir <- exp_descrip$src - hcst.dir <- exp_descrip$src - fcst.nmember <- exp_descrip$nmember$fcst - hcst.nmember <- exp_descrip$nmember$hcst - - ## TODO: it is necessary? - ##if ("accum" %in% names(reference_descrip)) { - ## accum <- unlist(reference_descrip$accum[store.freq][[1]]) - ##} else { - ## accum <- FALSE - ##} - - var_dir_obs <- reference_descrip[[store.freq]][variable] - var_dir_exp <- exp_descrip[[store.freq]][variable] - - # ----------- - obs.path <- paste0(archive$src, - obs.dir, store.freq, "/$var$", "$var_dir$", - "/$var$_$file_date$.nc") - - hcst.path <- paste0(archive$src, - hcst.dir, store.freq, "/$var$", "$var_dir$", - "$var$_$file_date$.nc") - - fcst.path <- paste0(archive$src, - hcst.dir, store.freq, "/$var$", "$var_dir$", - "/$var$_$file_date$.nc") - - # Define regrid parameters: - #------------------------------------------------------------------- - regrid_params <- get_regrid_params(recipe, archive) - - # Longitude circular sort and latitude check - #------------------------------------------------------------------- - circularsort <- check_latlon(lats.min, lats.max, lons.min, lons.max) - - if (recipe$Analysis$Variables$freq == "monthly_mean"){ - split_multiselected_dims = TRUE - } else { - split_multiselected_dims = FALSE - } - - # Load hindcast - #------------------------------------------------------------------- - hcst <- Start(dat = hcst.path, - var = variable, - var_dir = var_dir_exp, - file_date = sdates$hcst, - time = idxs$hcst, - var_dir_depends = 'var', - latitude = values(list(lats.min, lats.max)), - latitude_reorder = Sort(), - longitude = values(list(lons.min, lons.max)), - longitude_reorder = circularsort, - transform = regrid_params$fcst.transform, - transform_params = list(grid = regrid_params$fcst.gridtype, - method = regrid_params$fcst.gridmethod), - transform_vars = c('latitude', 'longitude'), - synonims = list(latitude = c('lat', 'latitude'), - longitude = c('lon', 'longitude'), - ensemble = c('member', 'ensemble')), - ensemble = indices(1:hcst.nmember), - metadata_dims = 'var', # change to just 'var'? - return_vars = list(latitude = 'dat', - longitude = 'dat', - time = 'file_date'), - split_multiselected_dims = split_multiselected_dims, - retrieve = TRUE) - - # Remove var_dir dimension - if ("var_dir" %in% names(dim(hcst))) { - hcst <- Subset(hcst, along = "var_dir", indices = 1, drop = "selected") - } - - if (store.freq %in% c("daily_mean", "daily")) { - # Adjusts dims for daily case, could be removed if startR allows - # multidim split - names(dim(hcst))[which(names(dim(hcst)) == 'file_date')] <- "syear" - default_dims <- c(dat = 1, var = 1, sday = 1, - sweek = 1, syear = 1, time = 1, - latitude = 1, longitude = 1, ensemble = 1) - default_dims[names(dim(hcst))] <- dim(hcst) - dim(hcst) <- default_dims - # Change time attribute dimensions - default_time_dims <- c(sday = 1, sweek = 1, syear = 1, time = 1) - names(dim(attr(hcst, "Variables")$common$time))[which(names( - dim(attr(hcst, "Variables")$common$time)) == 'file_date')] <- "syear" - default_time_dims[names(dim(attr(hcst, "Variables")$common$time))] <- - dim(attr(hcst, "Variables")$common$time) - dim(attr(hcst, "Variables")$common$time) <- default_time_dims - } - - # Convert hcst to s2dv_cube object - ## TODO: Give correct dimensions to $Dates - ## (sday, sweek, syear instead of file_date) - hcst <- as.s2dv_cube(hcst) - # Adjust dates for models where the time stamp goes into the next month - if (recipe$Analysis$Variables$freq == "monthly_mean") { - hcst$attrs$Dates[] <- hcst$attrs$Dates - seconds(exp_descrip$time_stamp_lag) - } - - # Load forecast - #------------------------------------------------------------------- - if (!is.null(recipe$Analysis$Time$fcst_year)) { - # the call uses file_date instead of fcst_syear so that it can work - # with the daily case and the current version of startR not allowing - # multiple dims split - - fcst <- Start(dat = fcst.path, - var = variable, - var_dir = var_dir_exp, - var_dir_depends = 'var', - file_date = sdates$fcst, - time = idxs$fcst, - latitude = values(list(lats.min, lats.max)), - latitude_reorder = Sort(), - longitude = values(list(lons.min, lons.max)), - longitude_reorder = circularsort, - transform = regrid_params$fcst.transform, - transform_params = list(grid = regrid_params$fcst.gridtype, - method = regrid_params$fcst.gridmethod), - transform_vars = c('latitude', 'longitude'), - synonims = list(latitude = c('lat', 'latitude'), - longitude = c('lon', 'longitude'), - ensemble = c('member', 'ensemble')), - ensemble = indices(1:fcst.nmember), - metadata_dims = 'var', - return_vars = list(latitude = 'dat', - longitude = 'dat', - time = 'file_date'), - split_multiselected_dims = split_multiselected_dims, - retrieve = TRUE) - - if ("var_dir" %in% names(dim(fcst))) { - fcst <- Subset(fcst, along = "var_dir", indices = 1, drop = "selected") - } - - if (store.freq %in% c("daily_mean", "daily")) { - # Adjusts dims for daily case, could be removed if startR allows - # multidim split - names(dim(fcst))[which(names(dim(fcst)) == 'file_date')] <- "syear" - default_dims <- c(dat = 1, var = 1, sday = 1, - sweek = 1, syear = 1, time = 1, - latitude = 1, longitude = 1, ensemble = 1) - default_dims[names(dim(fcst))] <- dim(fcst) - dim(fcst) <- default_dims - # Change time attribute dimensions - default_time_dims <- c(sday = 1, sweek = 1, syear = 1, time = 1) - names(dim(attr(fcst, "Variables")$common$time))[which(names( - dim(attr(fcst, "Variables")$common$time)) == 'file_date')] <- "syear" - default_time_dims[names(dim(attr(fcst, "Variables")$common$time))] <- - dim(attr(fcst, "Variables")$common$time) - dim(attr(fcst, "Variables")$common$time) <- default_time_dims - } - - # Convert fcst to s2dv_cube - fcst <- as.s2dv_cube(fcst) - # Adjust dates for models where the time stamp goes into the next month - if (recipe$Analysis$Variables$freq == "monthly_mean") { - fcst$attrs$Dates[] <- - fcst$attrs$Dates - seconds(exp_descrip$time_stamp_lag) - } - +Loading <- function(recipe) { + # Source correct function depending on time horizon + time_horizon <- tolower(recipe$Analysis$Horizon) + if (time_horizon == "seasonal") { + source("modules/Loading/R/load_seasonal.R") + data <- load_seasonal(recipe) + } else if (time_horizon == "decadal") { + source("modules/Loading/R/load_decadal.R") + data <- load_decadal(recipe) } else { - fcst <- NULL - } - - # Load reference - #------------------------------------------------------------------- - - # Obtain dates and date dimensions from the loaded hcst data to make sure - # the corresponding observations are loaded correctly. - dates <- hcst$attrs$Dates - dim(dates) <- hcst$dims[c("sday", "sweek", "syear", "time")] - - # Separate Start() call for monthly vs daily data - if (store.freq == "monthly_mean") { - - dates_file <- format(as.Date(dates, '%Y%m%d'), "%Y%m") - dim(dates_file) <- dim(dates) - - obs <- Start(dat = obs.path, - var = variable, - var_dir = var_dir_obs, - var_dir_depends = 'var', - file_date = dates_file, - latitude = values(list(lats.min, lats.max)), - latitude_reorder = Sort(), - longitude = values(list(lons.min, lons.max)), - longitude_reorder = circularsort, - transform = regrid_params$obs.transform, - transform_params = list(grid = regrid_params$obs.gridtype, - method = regrid_params$obs.gridmethod), - transform_vars = c('latitude', 'longitude'), - synonims = list(latitude = c('lat','latitude'), - longitude = c('lon','longitude')), - metadata_dims = 'var', - return_vars = list(latitude = 'dat', - longitude = 'dat', - time = 'file_date'), - split_multiselected_dims = TRUE, - retrieve = TRUE) - - } else if (store.freq %in% c("daily_mean", "daily")) { - - # Get year and month for file_date - dates_file <- sapply(dates, format, '%Y%m') - dim(dates_file) <- dim(dates) - # Set hour to 12:00 to ensure correct date retrieval for daily data - lubridate::hour(dates) <- 12 - lubridate::minute(dates) <- 00 - # Restore correct dimensions - dim(dates) <- dim(dates_file) - - obs <- Start(dat = obs.path, - var = variable, - var_dir = var_dir_obs, - var_dir_depends = 'var', - file_date = sort(unique(dates_file)), - time = dates, - time_var = 'time', - time_across = 'file_date', - merge_across_dims = TRUE, - merge_across_dims_narm = TRUE, - latitude = values(list(lats.min, lats.max)), - latitude_reorder = Sort(), - longitude = values(list(lons.min, lons.max)), - longitude_reorder = circularsort, - transform = regrid_params$obs.transform, - transform_params = list(grid = regrid_params$obs.gridtype, - method = regrid_params$obs.gridmethod), - transform_vars = c('latitude', 'longitude'), - synonims = list(latitude = c('lat','latitude'), - longitude = c('lon','longitude')), - metadata_dims = 'var', - return_vars = list(latitude = 'dat', - longitude = 'dat', - time = 'file_date'), - split_multiselected_dims = TRUE, - retrieve = TRUE) - } - - - # Remove var_dir dimension - if ("var_dir" %in% names(dim(obs))) { - obs <- Subset(obs, along = "var_dir", indices = 1, drop = "selected") - } - # Adds ensemble dim to obs (for consistency with hcst/fcst) - default_dims <- c(dat = 1, var = 1, sday = 1, - sweek = 1, syear = 1, time = 1, - latitude = 1, longitude = 1, ensemble = 1) - default_dims[names(dim(obs))] <- dim(obs) - dim(obs) <- default_dims - - # Convert obs to s2dv_cube - obs <- as.s2dv_cube(obs) - - # Check for consistency between hcst and obs grid - if (!(recipe$Analysis$Regrid$type == 'none')) { - if (!isTRUE(all.equal(as.vector(hcst$lat), as.vector(obs$lat)))) { - lat_error_msg <- paste("Latitude mismatch between hcst and obs.", - "Please check the original grids and the", - "regrid parameters in your recipe.") - error(recipe$Run$logger, lat_error_msg) - hcst_lat_msg <- paste0("First hcst lat: ", hcst$lat[1], - "; Last hcst lat: ", hcst$lat[length(hcst$lat)]) - info(recipe$Run$logger, hcst_lat_msg) - obs_lat_msg <- paste0("First obs lat: ", obs$lat[1], - "; Last obs lat: ", obs$lat[length(obs$lat)]) - info(recipe$Run$logger, obs_lat_msg) - stop("hcst and obs don't share the same latitudes.") - } - if (!isTRUE(all.equal(as.vector(hcst$lon), as.vector(obs$lon)))) { - lon_error_msg <- paste("Longitude mismatch between hcst and obs.", - "Please check the original grids and the", - "regrid parameters in your recipe.") - error(recipe$Run$logger, lon_error_msg) - hcst_lon_msg <- paste0("First hcst lon: ", hcst$lon[1], - "; Last hcst lon: ", hcst$lon[length(hcst$lon)]) - info(recipe$Run$logger, hcst_lon_msg) - obs_lon_msg <- paste0("First obs lon: ", obs$lon[1], - "; Last obs lon: ", obs$lon[length(obs$lon)]) - info(recipe$Run$logger, obs_lon_msg) - stop("hcst and obs don't share the same longitudes.") - - } + stop("Incorrect time horizon.") } - - # Remove negative values in accumulative variables - dictionary <- read_yaml("conf/variable-dictionary.yml") - for (var_idx in 1:length(variable)) { - var_name <- variable[var_idx] - if (dictionary$vars[[var_name]]$accum) { - info(recipe$Run$logger, - paste0("Accumulated variable ", var_name, - ": setting negative values to zero.")) - # obs$data[, var_idx, , , , , , , ] <- pmax(Subset(obs$data, - # along = "var", - # indices = var_idx, F), 0) - obs$data[, var_idx, , , , , , , ][obs$data[, var_idx, , , , , , , ] < 0] <- 0 - hcst$data[, var_idx, , , , , , , ][hcst$data[, var_idx, , , , , , , ] < 0] <- 0 - if (!is.null(fcst)) { - fcst$data[, var_idx, , , , , , , ][fcst$data[, var_idx, , , , , , , ] < 0] <- 0 - } - } - - # Convert prlr from m/s to mm/day - ## TODO: Make a unit conversion function - if (variable[[var_idx]] == "prlr") { - # Verify that the units are m/s and the same in obs and hcst - if (((obs$attrs$Variable$metadata[[var_name]]$units == "m s-1") || - (obs$attrs$Variable$metadata[[var_name]]$units == "m s**-1")) && - ((hcst$attrs$Variable$metadata[[var_name]]$units == "m s-1") || - (hcst$attrs$Variable$metadata[[var_name]]$units == "m s**-1"))) { - info(recipe$Run$logger, "Converting precipitation from m/s to mm/day.") - obs$data[, var_idx, , , , , , , ] <- - obs$data[, var_idx, , , , , , , ]*86400*1000 - obs$attrs$Variable$metadata[[var_name]]$units <- "mm/day" - hcst$data[, var_idx, , , , , , , ] <- - hcst$data[, var_idx, , , , , , , ]*86400*1000 - hcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day" - if (!is.null(fcst)) { - fcst$data[, var_idx, , , , , , , ] <- - fcst$data[, var_idx, , , , , , , ]*86400*1000 - fcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day" - } - } - } - } - # Compute anomalies if requested - # Print a summary of the loaded data for the user, for each object - if (recipe$Run$logger$threshold <= 2) { - data_summary(hcst, recipe) - data_summary(obs, recipe) - if (!is.null(fcst)) { - data_summary(fcst, recipe) - } - } - - info(recipe$Run$logger, - "##### DATA LOADING COMPLETED SUCCESSFULLY #####") - - ############################################################################ - # - # CHECKS ON MISSING FILES - # - ############################################################################ - - #obs.NA_dates.ind <- Apply(obs, - # fun=(function(x){ all(is.na(x))}), - # target_dims=c('time', 'latitude', 'longitude'))[[1]] - #obs.NA_dates <- dates_file[obs.NA_dates.ind] - #obs.NA_dates <- obs.NA_dates[order(obs.NA_dates)] - #obs.NA_files <- paste0(obs.dir, store.freq,"/",variable,"_", - # freq.obs,"obs.grid","/",variable,"_",obs.NA_dates,".nc") - # - #if (any(is.na(hcst))){ - # fatal(recipe$Run$logger, - # paste(" ERROR: MISSING HCST VALUES FOUND DURING LOADING # ", - # " ################################################# ", - # " ###### MISSING FILES #### ", - # " ################################################# ", - # "hcst files:", - # hcst.NA_files, - # " ################################################# ", - # " ################################################# ", - # sep="\n")) - # quit(status = 1) - #} - # - #if (any(is.na(obs)) && !identical(obs.NA_dates,character(0))){ - # fatal(recipe$logger, - # paste(" ERROR: MISSING OBS VALUES FOUND DURING LOADING # ", - # " ################################################# ", - # " ###### MISSING FILES #### ", - # " ################################################# ", - # "obs files:", - # obs.NA_files, - # " ################################################# ", - # " ################################################# ", - # sep="\n")) - # quit(status=1) - #} - # - #info(recipe$logger, - # "######### DATA LOADING COMPLETED SUCCESFULLY ##############") - - ############################################################################ - ############################################################################ - - return(list(hcst = hcst, fcst = fcst, obs = obs)) - + return(data) } + diff --git a/modules/Saving/Saving.R b/modules/Saving/Saving.R index 2d067ba3..7d0520c2 100644 --- a/modules/Saving/Saving.R +++ b/modules/Saving/Saving.R @@ -14,11 +14,12 @@ source("modules/Saving/R/get_time.R") source("modules/Saving/R/get_latlon.R") source("modules/Saving/R/get_global_attributes.R") source("modules/Saving/R/drop_dims.R") -save_data <- function(recipe, data, - skill_metrics = NULL, - probabilities = NULL, - agg = 'global', - archive = NULL) { + +Saving <- function(recipe, data, + skill_metrics = NULL, + probabilities = NULL, + agg = 'global', + archive = NULL) { # Wrapper for the saving functions. # recipe: The auto-s2s recipe # archive: The auto-s2s archive diff --git a/modules/Skill/Skill.R b/modules/Skill/Skill.R index ecb89224..c6ad37bd 100644 --- a/modules/Skill/Skill.R +++ b/modules/Skill/Skill.R @@ -57,7 +57,7 @@ source("modules/Skill/R/CRPS_clim.R") # compute_skill_metrics <- function(recipe, data$hcst, obs, # clim_data$hcst = NULL, # clim_obs = NULL) { -compute_skill_metrics <- function(recipe, data, agg = 'global') { +Skill <- function(recipe, data, agg = 'global') { # data$hcst: s2dv_cube containing the hindcast # obs: s2dv_cube containing the observations @@ -343,7 +343,7 @@ compute_skill_metrics <- function(recipe, data, agg = 'global') { return(skill_metrics) } -compute_probabilities <- function(recipe, data) { +Probabilities <- function(recipe, data) { ## TODO: Do hcst and fcst at the same time if (is.null(recipe$Analysis$ncores)) { diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R index 03a4d5a2..228b01c2 100644 --- a/modules/Visualization/Visualization.R +++ b/modules/Visualization/Visualization.R @@ -11,11 +11,11 @@ source("modules/Visualization/R/tmp/PlotRobinson.R") source("modules/Visualization/R/plot_most_likely_terciles_map.R") source("modules/Visualization/R/plot_ensemble_mean.R") -plot_data <- function(recipe, - data, - skill_metrics = NULL, - probabilities = NULL, - significance = F) { +Visualization <- function(recipe, + data, + skill_metrics = NULL, + probabilities = NULL, + significance = F) { # Try to produce and save several basic plots. # recipe: the auto-s2s recipe as read by read_yaml() # data: list containing the hcst, obs and (optional) fcst s2dv_cube objects -- GitLab From f3c702a5f69c85aa712e8f21db5cf1dcde9bc742 Mon Sep 17 00:00:00 2001 From: Victoria Agudetse Roures Date: Wed, 5 Jul 2023 16:56:19 +0200 Subject: [PATCH 3/3] Change function names in unit tests and source libs.R in Loading.R --- modules/Loading/Loading.R | 2 ++ modules/Loading/R/load_decadal.R | 1 - modules/Loading/R/load_seasonal.R | 2 -- modules/test_parallel_workflow.R | 16 ++++++++-------- tests/testthat/test-decadal_daily_1.R | 10 +++++----- tests/testthat/test-decadal_monthly_1.R | 24 ++++++++++++------------ tests/testthat/test-decadal_monthly_2.R | 16 ++++++++-------- tests/testthat/test-decadal_monthly_3.R | 10 +++++----- tests/testthat/test-seasonal_NAO.R | 6 +++--- tests/testthat/test-seasonal_daily.R | 6 +++--- tests/testthat/test-seasonal_monthly.R | 18 +++++++++--------- 11 files changed, 55 insertions(+), 56 deletions(-) diff --git a/modules/Loading/Loading.R b/modules/Loading/Loading.R index 9b9ee86b..78661bf2 100644 --- a/modules/Loading/Loading.R +++ b/modules/Loading/Loading.R @@ -1,3 +1,5 @@ +source("tools/libs.R") + Loading <- function(recipe) { # Source correct function depending on time horizon time_horizon <- tolower(recipe$Analysis$Horizon) diff --git a/modules/Loading/R/load_decadal.R b/modules/Loading/R/load_decadal.R index 7d0b2d8b..6db430d7 100644 --- a/modules/Loading/R/load_decadal.R +++ b/modules/Loading/R/load_decadal.R @@ -11,7 +11,6 @@ source("modules/Loading/helper_loading_decadal.R") source("modules/Loading/R/dates2load.R") source("modules/Loading/R/check_latlon.R") source("modules/Loading/R/get_timeidx.R") -source("tools/libs.R") #==================================================================== diff --git a/modules/Loading/R/load_seasonal.R b/modules/Loading/R/load_seasonal.R index c54ef908..27f654ab 100644 --- a/modules/Loading/R/load_seasonal.R +++ b/modules/Loading/R/load_seasonal.R @@ -4,8 +4,6 @@ source("/esarchive/scratch/vagudets/repos/csoperational/R/get_regrid_params.R") source("modules/Loading/R/dates2load.R") source("modules/Loading/R/get_timeidx.R") source("modules/Loading/R/check_latlon.R") -## TODO: Move to prepare_outputs.R -source("tools/libs.R") load_seasonal <- function(recipe) { diff --git a/modules/test_parallel_workflow.R b/modules/test_parallel_workflow.R index 30c9cb91..28da1f30 100644 --- a/modules/test_parallel_workflow.R +++ b/modules/test_parallel_workflow.R @@ -9,17 +9,17 @@ args = commandArgs(trailingOnly = TRUE) recipe_file <- args[1] recipe <- read_atomic_recipe(recipe_file) # Load datasets -data <- load_datasets(recipe) +data <- Loading(recipe) # Calibrate datasets -data <- calibrate_datasets(recipe, data) +data <- Calibration(recipe, data) # Compute anomalies -data <- compute_anomalies(recipe, data) +data <- Anomalies(recipe, data) # Compute skill metrics -skill_metrics <- compute_skill_metrics(recipe, data) +skill_metrics <- Skill(recipe, data) # Compute percentiles and probability bins -probabilities <- compute_probabilities(recipe, data) +probabilities <- Probabilities(recipe, data) # Export all data to netCDF -save_data(recipe, data, skill_metrics, probabilities) +Saving(recipe, data, skill_metrics, probabilities) # Plot data -plot_data(recipe, data, skill_metrics, probabilities, - significance = T) +Visualization(recipe, data, skill_metrics, probabilities, + significance = T) diff --git a/tests/testthat/test-decadal_daily_1.R b/tests/testthat/test-decadal_daily_1.R index c26b8978..c847fd10 100644 --- a/tests/testthat/test-decadal_daily_1.R +++ b/tests/testthat/test-decadal_daily_1.R @@ -2,7 +2,7 @@ context("Decadal daily data - 1") ########################################### -source("modules/Loading/Loading_decadal.R") +source("modules/Loading/Loading.R") source("modules/Calibration/Calibration.R") source("modules/Skill/Skill.R") source("modules/Saving/Saving.R") @@ -13,18 +13,18 @@ archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive_decadal.yml"))$ar # Load datasets suppressWarnings({invisible(capture.output( -data <- load_datasets(recipe) +data <- Loading(recipe) ))}) ## Calibrate datasets #suppressWarnings({invisible(capture.output( -# calibrated_data <- calibrate_datasets(data, recipe) +# calibrated_data <- Calibration(data, recipe) #))}) # ## Compute skill metrics #suppressWarnings({invisible(capture.output( -#skill_metrics <- compute_skill_metrics(calibrated_data$hcst, data$obs, -# recipe, na.rm = T, ncores = 4) +#skill_metrics <- Skill(calibrated_data$hcst, data$obs, +# recipe, na.rm = T, ncores = 4) #))}) #====================================== diff --git a/tests/testthat/test-decadal_monthly_1.R b/tests/testthat/test-decadal_monthly_1.R index a5c95e0c..3346e529 100644 --- a/tests/testthat/test-decadal_monthly_1.R +++ b/tests/testthat/test-decadal_monthly_1.R @@ -2,7 +2,7 @@ context("Decadal monthly data - 1") ########################################### -source("modules/Loading/Loading_decadal.R") +source("modules/Loading/Loading.R") source("modules/Calibration/Calibration.R") source("modules/Skill/Skill.R") source("modules/Saving/Saving.R") @@ -14,27 +14,27 @@ archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive_decadal.yml"))$ar # Load datasets suppressWarnings({invisible(capture.output( -data <- load_datasets(recipe) +data <- Loading(recipe) ))}) # Calibrate datasets suppressWarnings({invisible(capture.output( - calibrated_data <- calibrate_datasets(recipe, data) + calibrated_data <- Calibration(recipe, data) ))}) # Compute skill metrics suppressWarnings({invisible(capture.output( -skill_metrics <- compute_skill_metrics(recipe, calibrated_data) +skill_metrics <- Skill(recipe, calibrated_data) ))}) suppressWarnings({invisible(capture.output( -probs <- compute_probabilities(recipe, calibrated_data) +probs <- Probabilities(recipe, calibrated_data) ))}) # Plotting suppressWarnings({invisible(capture.output( -plot_data(recipe = recipe, data = calibrated_data, - skill_metrics = skill_metrics, probabilities = probs, - significance = T) +Visualization(recipe = recipe, data = calibrated_data, + skill_metrics = skill_metrics, probabilities = probs, + significance = T) ))}) #====================================== @@ -297,20 +297,20 @@ archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive_decadal.yml"))$ar # Load datasets suppressWarnings({invisible(capture.output( -data_b <- load_datasets(recipe) +data_b <- Loading(recipe) ))}) # Calibrate datasets suppressWarnings({invisible(capture.output( - calibrated_data_b <- calibrate_datasets(recipe, data_b) + calibrated_data_b <- Calibration(recipe, data_b) ))}) # Compute skill metrics suppressWarnings({invisible(capture.output( -skill_metrics_b <- compute_skill_metrics(recipe, calibrated_data_b) +skill_metrics_b <- Skill(recipe, calibrated_data_b) ))}) suppressWarnings({invisible(capture.output( -probs_b <- compute_probabilities(recipe, calibrated_data_b) +probs_b <- Probabilities(recipe, calibrated_data_b) ))}) diff --git a/tests/testthat/test-decadal_monthly_2.R b/tests/testthat/test-decadal_monthly_2.R index 2b605109..9adc16b6 100644 --- a/tests/testthat/test-decadal_monthly_2.R +++ b/tests/testthat/test-decadal_monthly_2.R @@ -2,7 +2,7 @@ context("Decadal monthly data - 2") ########################################### -source("modules/Loading/Loading_decadal.R") +source("modules/Loading/Loading.R") source("modules/Calibration/Calibration.R") source("modules/Skill/Skill.R") source("modules/Saving/Saving.R") @@ -13,27 +13,27 @@ recipe <- prepare_outputs(recipe_file) # Load datasets suppressWarnings({invisible(capture.output( -data <- load_datasets(recipe) +data <- Loading(recipe) ))}) # Calibrate datasets suppressWarnings({invisible(capture.output( - calibrated_data <- calibrate_datasets(recipe, data) + calibrated_data <- Calibration(recipe, data) ))}) # Compute skill metrics suppressMessages({invisible(capture.output( -skill_metrics <- compute_skill_metrics(recipe, calibrated_data) +skill_metrics <- Skill(recipe, calibrated_data) ))}) suppressWarnings({invisible(capture.output( -probs <- compute_probabilities(recipe, calibrated_data) +probs <- Probabilities(recipe, calibrated_data) ))}) # Plotting suppressWarnings({invisible(capture.output( -plot_data(recipe = recipe, data = calibrated_data, - skill_metrics = skill_metrics, - probabilities = probs, significance = T) +Visualization(recipe = recipe, data = calibrated_data, + skill_metrics = skill_metrics, + probabilities = probs, significance = T) ))}) diff --git a/tests/testthat/test-decadal_monthly_3.R b/tests/testthat/test-decadal_monthly_3.R index 7232ebfd..d1b42fd5 100644 --- a/tests/testthat/test-decadal_monthly_3.R +++ b/tests/testthat/test-decadal_monthly_3.R @@ -2,7 +2,7 @@ context("Decadal monthly data - 3") ########################################### -source("modules/Loading/Loading_decadal.R") +source("modules/Loading/Loading.R") source("modules/Calibration/Calibration.R") source("modules/Skill/Skill.R") source("modules/Saving/Saving.R") @@ -13,20 +13,20 @@ archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive_decadal.yml"))$ar # Load datasets suppressWarnings({invisible(capture.output( -data <- load_datasets(recipe) +data <- Loading(recipe) ))}) # Calibrate datasets suppressWarnings({invisible(capture.output( - calibrated_data <- calibrate_datasets(recipe, data) + calibrated_data <- Calibration(recipe, data) ))}) # Compute skill metrics suppressWarnings({invisible(capture.output( -skill_metrics <- compute_skill_metrics(recipe, calibrated_data) +skill_metrics <- Skill(recipe, calibrated_data) ))}) suppressWarnings({invisible(capture.output( -probs <- compute_probabilities(recipe, calibrated_data) +probs <- Probabilities(recipe, calibrated_data) ))}) #====================================== diff --git a/tests/testthat/test-seasonal_NAO.R b/tests/testthat/test-seasonal_NAO.R index 7a03ecbf..7a9c9de8 100644 --- a/tests/testthat/test-seasonal_NAO.R +++ b/tests/testthat/test-seasonal_NAO.R @@ -12,12 +12,12 @@ archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive.yml"))$archive # Load datasets suppressWarnings({invisible(capture.output( -data <- load_datasets(recipe) +data <- Loading(recipe) ))}) # Compute anomalies suppressWarnings({invisible(capture.output( -ano_data <- compute_anomalies(recipe, data) +ano_data <- Anomalies(recipe, data) ))}) # Compute index NAO @@ -27,7 +27,7 @@ nao_data <- Indices(recipe = recipe, data = ano_data) # Compute skill metrics suppressWarnings({invisible(capture.output( -skill_metrics <- compute_skill_metrics(recipe, nao_data, agg = 'region') +skill_metrics <- Skill(recipe, nao_data, agg = 'region') ))}) diff --git a/tests/testthat/test-seasonal_daily.R b/tests/testthat/test-seasonal_daily.R index 8a51cc86..6cfa4384 100644 --- a/tests/testthat/test-seasonal_daily.R +++ b/tests/testthat/test-seasonal_daily.R @@ -9,17 +9,17 @@ recipe_file <- "tests/recipes/recipe-seasonal_daily_1.yml" recipe <- prepare_outputs(recipe_file, disable_checks = F) # Load datasets suppressWarnings({invisible(capture.output( -data <- load_datasets(recipe) +data <- Loading(recipe) ))}) # Calibrate data suppressWarnings({invisible(capture.output( -calibrated_data <- calibrate_datasets(recipe, data) +calibrated_data <- Calibration(recipe, data) ))}) # Compute skill metrics suppressWarnings({invisible(capture.output( -skill_metrics <- compute_skill_metrics(recipe, calibrated_data) +skill_metrics <- Skill(recipe, calibrated_data) ))}) test_that("1. Loading", { diff --git a/tests/testthat/test-seasonal_monthly.R b/tests/testthat/test-seasonal_monthly.R index 0e311424..d112acf9 100644 --- a/tests/testthat/test-seasonal_monthly.R +++ b/tests/testthat/test-seasonal_monthly.R @@ -12,34 +12,34 @@ archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive.yml"))$archive # Load datasets suppressWarnings({invisible(capture.output( -data <- load_datasets(recipe) +data <- Loading(recipe) ))}) # Calibrate data suppressWarnings({invisible(capture.output( -calibrated_data <- calibrate_datasets(recipe, data) +calibrated_data <- Calibration(recipe, data) ))}) # Compute skill metrics suppressWarnings({invisible(capture.output( -skill_metrics <- compute_skill_metrics(recipe, calibrated_data) +skill_metrics <- Skill(recipe, calibrated_data) ))}) suppressWarnings({invisible(capture.output( -probs <- compute_probabilities(recipe, calibrated_data) +probs <- Probabilities(recipe, calibrated_data) ))}) # Saving suppressWarnings({invisible(capture.output( -save_data(recipe = recipe, data = calibrated_data, - skill_metrics = skill_metrics, probabilities = probs) +Saving(recipe = recipe, data = calibrated_data, + skill_metrics = skill_metrics, probabilities = probs) ))}) # Plotting suppressWarnings({invisible(capture.output( -plot_data(recipe = recipe, data = calibrated_data, - skill_metrics = skill_metrics, probabilities = probs, - significance = T) +Visualization(recipe = recipe, data = calibrated_data, + skill_metrics = skill_metrics, probabilities = probs, + significance = T) ))}) outdir <- get_dir(recipe = recipe, variable = data$hcst$attrs$Variable$varName) -- GitLab