diff --git a/conf/archive.yml b/conf/archive.yml index 9d994ee511af7fe1646f355806a8bade6026d107..c50909c3ed556d255c81c23370a0322e90694f14 100644 --- a/conf/archive.yml +++ b/conf/archive.yml @@ -19,7 +19,7 @@ archive: calendar: "proleptic_gregorian" reference_grid: "/esarchive/exp/ecmwf/system5c3s/monthly_mean/tas_f6h/tas_20180501.nc" system7c3s: - name: "Méteo-France System 7" + name: "Meteo-France System 7" institution: "European Centre for Medium-Range Weather Forecasts" src: "exp/meteofrance/system7c3s/" monthly_mean: {"tas":"_f6h/", "g500":"_f12h/", diff --git a/modules/Calibration/Calibration.R b/modules/Calibration/Calibration.R index 85b1b0070612a4e96f9cf75652c2f56954491a00..15e51e82ed02218abad394ecc678126799fc1d0c 100644 --- a/modules/Calibration/Calibration.R +++ b/modules/Calibration/Calibration.R @@ -1,7 +1,7 @@ ## TODO: Remove once Alba's fun is merged in CSTools source("tools/tmp/CST_Calibration.R") - +source("tools/tmp/CST_QuantileMapping.R") ## Entry params data and recipe? calibrate_datasets <- function(data, recipe) { # Function that calibrates the hindcast using the method stated in the @@ -93,29 +93,32 @@ calibrate_datasets <- function(data, recipe) { "frequency. Only quantile mapping 'qmap' is implemented.") } # Calibrate the hindcast + dim_order <- names(dim(data$hcst$data)) hcst_calibrated <- CST_QuantileMapping(data$hcst, data$obs, exp_cor = NULL, - sample_dims = c("syear", - "time", - "ensemble"), - sample_length = NULL, + sdate_dim = "syear", + memb_dim = "ensemble", + # window_dim = "time", method = "QUANT", - wet.day = FALSE, ncores = ncores, - na.rm = na.rm) + na.rm = na.rm, + wet.day = F) + # Restore dimension order + hcst_calibrated$data <- Reorder(hcst_calibrated$data, dim_order) if (!is.null(data$fcst)) { # Calibrate the forecast fcst_calibrated <- CST_QuantileMapping(data$hcst, data$obs, exp_cor = data$fcst, - sample_dims = c("syear", - "time", - "ensemble"), - sample_length = NULL, + sdate_dim = "syear", + memb_dim = "ensemble", + # window_dim = "time", method = "QUANT", - wet.day = FALSE, ncores = ncores, - na.rm = na.rm) + na.rm = na.rm, + wet.day = F) + # Restore dimension order + fcst_calibrated$data <- Reorder(fcst_calibrated$data, dim_order) } else { fcst_calibrated <- NULL } diff --git a/tests/testthat/test-seasonal_daily.R b/tests/testthat/test-seasonal_daily.R index c37b551486bb65a409079e762bb5145740ed4f0a..237674e09715273e3890961c3dd942103f996d3b 100644 --- a/tests/testthat/test-seasonal_daily.R +++ b/tests/testthat/test-seasonal_daily.R @@ -123,17 +123,17 @@ c(dat = 1, var = 1, sday = 1, sweek = 1, syear = 4, time = 31, latitude = 4, lon ) expect_equal( mean(calibrated_data$hcst$data), -289.6468, +289.5354, tolerance = 0.0001 ) expect_equal( as.vector(drop(calibrated_data$hcst$data)[1, 1:4, 2, 3, 4]), -c(295.1077, 294.2161, 294.5801, 292.6326), +c(291.5555, 291.9029, 293.2685, 290.7782), tolerance = 0.0001 ) expect_equal( range(calibrated_data$hcst$data), -c(283.9447, 297.7496), +c(284.2823, 296.7545), tolerance = 0.0001 ) }) @@ -160,7 +160,7 @@ c(time = 31, latitude = 4, longitude = 4) ) expect_equal( skill_metrics$enscorr_specs[1:3, 1, 1], -c(0.8159317, 0.8956195, 0.8355627), +c(0.7509920, 0.6514916, 0.5118371), tolerance=0.0001 ) }) diff --git a/tools/tmp/CST_QuantileMapping.R b/tools/tmp/CST_QuantileMapping.R new file mode 100644 index 0000000000000000000000000000000000000000..a4e4a9d601d504a27edc5128d571796ee45e87bf --- /dev/null +++ b/tools/tmp/CST_QuantileMapping.R @@ -0,0 +1,325 @@ +#'Quantile Mapping for seasonal or decadal forecast data +#' +#'@description This function is a wrapper of fitQmap and doQmap from package +#''qmap' to be applied on the object of class 's2dv_cube'. The quantile mapping +#'adjustment between an experiment, typically a hindcast, and observation is +#'applied to the experiment itself or to a provided forecast. + +#'@author Nuria Perez-Zanon, \email{nuria.perez@bsc.es} +#'@param exp An object of class \code{s2dv_cube}. +#'@param obs An object of class \code{s2dv_cube}. +#'@param exp_cor An object of class \code{s2dv_cube} in which the quantile +#' mapping correction should be applied. If it is not specified, the correction +#' is applied in object 'exp'. +#'@param sdate_dim A character string indicating the dimension name in which +#' cross-validation would be applied when exp_cor is not provided. 'sdate' by +#' default. +#'@param memb_dim A character string indicating the dimension name where +#' ensemble members are stored in the experimental arrays. 'member' by default. +#'@param window_dim A character string indicating the dimension name where +#' samples have been stored. It can be NULL (default) in case all samples are +#' used. +#'@param method A character string indicating the method to be used:'PTF', +#' 'DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping +#' 'QUANT' is used. +#'@param na.rm A logical value indicating if missing values should be removed +#' (FALSE by default). +#'@param ncores An integer indicating the number of cores for parallel +#' computation using multiApply function. The default value is NULL (1). +#'@param ... Additional parameters to be used by the method choosen. See qmap +#' package for details. +#' +#'@return An object of class \code{s2dv_cube} containing the experimental data +#'after applying the quantile mapping correction. +#' +#'@seealso \code{\link[qmap]{fitQmap}} and \code{\link[qmap]{doQmap}} +#'@examples +#'# Use synthetic data +#'exp <- NULL +#'exp$data <- 1 : c(1 * 10 * 20 * 60 * 6 * 7) +#'dim(exp$data) <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , +#' lat = 6, lon = 7) +#'class(exp) <- 's2dv_cube' +#'obs$data <- 101 : c(100 + 1 * 1 * 20 * 60 * 6 * 7) +#'dim(obs$data) <- c(dataset = 1, member = 1, sdate = 20, ftime = 60 , +#' lat = 6, lon = 7) +#'class(obs) <- 's2dv_cube' +#'res <- CST_QuantileMapping(exp, obs) +#' +#'# Use data in package +#'exp <- lonlat_temp$exp +#'exp$data <- exp$data[, , 1:4, , 1:2, 1:3] +#'dim(exp$data) <- c(dataset = 1, member = 15, sdate = 4, ftime = 3, +#' lat = 2, lon = 3) +#'obs <- lonlat_temp$obs +#'obs$data <- obs$data[, , 1:4, , 1:2, 1:3] +#'dim(obs$data) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3, +#' lat = 2, lon = 3) +#'exp_cor <- lonlat_temp$exp +#'exp_cor$data <- exp_cor$data[, 1, 5:6, , 1:2, 1:3] +#'dim(exp_cor$data) <- c(dataset = 1, member = 1, sdate = 2, ftime = 3, +#' lat = 2, lon = 3) +#'res <- CST_QuantileMapping(exp, obs, exp_cor, window_dim = 'ftime') +#' +#'@import qmap +#'@import multiApply +#'@import s2dv +#'@export +CST_QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', + memb_dim = 'member', window_dim = NULL, + method = 'QUANT', na.rm = FALSE, + ncores = NULL, ...) { + if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) { + stop("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + if (!is.null(exp_cor)) { + if (!inherits(exp_cor, 's2dv_cube')) { + stop("Parameter 'exp_cor' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + } + + QMapped <- QuantileMapping(exp = exp$data, obs = obs$data, + exp_cor = exp_cor$data, + sdate_dim = sdate_dim, memb_dim = memb_dim, + window_dim = window_dim, method = method, + na.rm = na.rm, ncores = ncores, ...) + if (is.null(exp_cor)) { + exp$data <- QMapped + exp$Datasets <- c(exp$Datasets, obs$Datasets) + exp$source_files <- c(exp$source_files, obs$source_files) + return(exp) + + } else { + exp_cor$data <- QMapped + exp_cor$Datasets <- c(exp_cor$Datasets, exp$Datasets, obs$Datasets) + exp_cor$source_files <- c(exp_cor$source_files, exp$source_files, obs$source_files) + return(exp_cor) + } + + +} + +#'Quantile Mapping for seasonal or decadal forecast data +#' +#'@description This function is a wrapper of fitQmap and doQmap from package +#''qmap' to be applied on multi-dimensional arrays. The quantile mapping +#'adjustment between an experiment, typically a hindcast, and observation is +#'applied to the experiment itself or to a provided forecast. +#' +#'@author Nuria Perez-Zanon, \email{nuria.perez@bsc.es} +#'@param exp A multidimensional array with named dimensions containing the +#' hindcast. +#'@param obs A multidimensional array with named dimensions containing the +#' reference dataset. +#'@param exp_cor A multidimensional array with named dimensions in which the +#' quantile mapping correction should be applied. If it is not specified, the +#' correction is applied on object 'exp'. +#'@param sdate_dim A character string indicating the dimension name in which +#' cross-validation would be applied when exp_cor is not provided. 'sdate' by +#' default. +#'@param memb_dim A character string indicating the dimension name where +#' ensemble members are stored in the experimental arrays. 'member' by +#' default. +#'@param window_dim A character string indicating the dimension name where +#' samples have been stored. It can be NULL (default) in case all samples are +#' used. +#'@param method A character string indicating the method to be used: 'PTF', +#' 'DIST','RQUANT','QUANT','SSPLIN'. By default, the empirical quantile mapping +#' 'QUANT' is used. +#'@param na.rm A logical value indicating if missing values should be removed +#' (FALSE by default). +#'@param ncores An integer indicating the number of cores for parallel +#' computation using multiApply function. The default value is NULL (1). +#'@param ... Additional parameters to be used by the method choosen. See qmap +#' package for details. +#' +#'@return An array containing the experimental data after applying the quantile +#'mapping correction. +#' +#'@seealso \code{\link[qmap]{fitQmap}} and \code{\link[qmap]{doQmap}} +#'@examples +#'# Use synthetic data +#'exp <- 1 : c(1 * 10 * 20 * 60 * 6 * 7) +#'dim(exp) <- c(dataset = 1, member = 10, sdate = 20, ftime = 60 , +#' lat = 6, lon = 7) +#'class(exp) <- 's2dv_cube' +#'obs <- 101 : c(100 + 1 * 1 * 20 * 60 * 6 * 7) +#'dim(obs) <- c(dataset = 1, member = 1, sdate = 20, ftime = 60 , +#' lat = 6, lon = 7) +#'res <- QuantileMapping(exp, obs) +#'# Use data in package +#'exp <- lonlat_temp$exp$data[, , 1:4, , 1:2, 1:3] +#'dim(exp) <- c(dataset = 1, member = 15, sdate = 4, ftime = 3, +#' lat = 2, lon = 3) +#'obs <- lonlat_temp$obs$data[, , 1:4, , 1:2, 1:3] +#'dim(obs) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3, +#' lat = 2, lon = 3) +#'exp_cor <- lonlat_temp$exp$data[, 1, 5:6, , 1:2, 1:3] +#'dim(exp_cor) <- c(dataset = 1, member = 1, sdate = 2, ftime = 3, +#' lat = 2, lon = 3) +#'res <- QuantileMapping(exp, obs, exp_cor, window_dim = 'ftime') +#' +#'@import qmap +#'@import multiApply +#'@import s2dv +#'@export +QuantileMapping <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', + memb_dim = 'member', window_dim = NULL, + method = 'QUANT', + na.rm = FALSE, ncores = NULL, ...) { + # exp and obs + obsdims <- names(dim(obs)) + expdims <- names(dim(exp)) + if (!is.array(exp) || !is.numeric(exp)) { + stop("Parameter 'exp' must be a numeric array.") + } + if (!is.array(obs) || !is.numeric(obs)) { + stop("Parameter 'obs' must be a numeric array.") + } + if (is.null(expdims)) { + stop("Parameter 'exp' must have dimension names.") + } + if (is.null(obsdims)) { + stop("Parameter 'obs' must have dimension names.") + } + # sdate_dim + if (!is.character(sdate_dim) | length(sdate_dim) != 1) { + stop("Parameter 'sdate_dim' must be a character string.") + } + if (!sdate_dim %in% expdims | !sdate_dim %in% obsdims) { + stop("Parameter 'sdate_dim' is not found in 'exp' or 'obs' dimension.") + } + if (dim(exp)[sdate_dim] == 1 || dim(obs)[sdate_dim] == 1) { + stop("Parameter 'exp' and 'obs' must have dimension length of 'sdate_dim' bigger than 1.") + } + # exp_cor + if (!is.null(exp_cor)) { + if (is.null(names(dim(exp_cor)))) { + stop("Parameter 'exp_cor' must have dimension names.") + } + if (!sdate_dim %in% names(dim(exp_cor))) { + stop("Parameter 'sdate_dim' is not found in 'exp_cor' dimension.") + } + } + # method + if (!(method %in% c('PTF', 'DIST', 'RQUANT', 'QUANT', 'SSPLIN')) | length(method) != 1) { + stop("Parameter 'method' must be one of the following methods: ", + "'PTF', 'DIST', 'RQUANT', 'QUANT', 'SSPLIN'.") + } + # memb_dim + if (!all(memb_dim %in% obsdims)) { + obs <- InsertDim(obs, posdim = 1, lendim = 1, + name = memb_dim[!(memb_dim %in% obsdims)]) + } + if (any(!memb_dim %in% expdims)) { + stop("Parameter 'memb_dim' is not found in 'exp' dimensions.") + } + sample_dims <- c(memb_dim, sdate_dim) + # window_dim + if (!is.null(window_dim)) { + if (!(window_dim %in% obsdims)) { + stop("Parameter 'window_dim' is not found in 'obs'.") + } + obs <- CSTools::MergeDims(obs, c(memb_dim, window_dim)) + if (window_dim %in% expdims) { + exp <- CSTools::MergeDims(exp, c(memb_dim, window_dim)) + warning("Parameter 'window_dim' is found in exp and is merged to 'memb_dim'.") + } + } + # na.rm + if (!is.logical(na.rm) | length(na.rm) > 1) { + stop("Parameter 'na.rm' must be one logical value.") + } + # ncores + if (!is.null(ncores)) { + if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 | + length(ncores) > 1) { + stop("Parameter 'ncores' must be either NULL or a positive integer.") + } + } + + ############################### + + if (!is.null(exp_cor)) { + qmaped <- Apply(list(exp, obs, exp_cor), target_dims = sample_dims, + fun = .qmapcor, method = method, sdate_dim = sdate_dim, + na.rm = na.rm, ..., + ncores = ncores)$output1 + } else { + qmaped <- Apply(list(exp, obs), target_dims = sample_dims, + fun = .qmapcor, exp_cor = NULL, method = method, + sdate_dim = sdate_dim, na.rm = na.rm, ..., + ncores = ncores)$output1 + } + return(qmaped) +} + +.qmapcor <- function(exp, obs, exp_cor = NULL, sdate_dim = 'sdate', method = 'QUANT', + na.rm = FALSE, ...) { + + # exp: [memb (+ window), sdate] + # obs: [memb (+ window), sdate] + # exp_cor: NULL or [memb, sdate] + + if (is.null(exp_cor)) { + applied <- exp * NA + for (sd in 1:dim(exp)[sdate_dim]) { + if (na.rm) { + # select start date for cross-val + nas_pos <- which(!is.na(exp[, sd])) + obs2 <- as.vector(obs[, -sd]) + exp2 <- as.vector(exp[, -sd]) + exp_cor2 <- as.vector(exp[, sd]) + # remove NAs + obs2 <- obs2[!is.na(obs2)] + exp2 <- exp2[!is.na(exp2)] + exp_cor2 <- exp_cor2[!is.na(exp_cor2)] + tryCatch({ + adjust <- fitQmap(obs2, exp2, method = method, ...) + applied[nas_pos, sd] <- doQmap(exp_cor2, adjust, ...) + }, + error = function(error_message) { + return(applied[, sd]) + }) + } else { + # na.rm = FALSE shouldn't fail, just return NA + if (anyNA(obs[, -sd]) | anyNA(exp[, -sd])) { + applied[, sd] <- NA + } else { + adjust <- fitQmap(as.vector(obs[, -sd]), as.vector(exp[, -sd]), + method = method, ...) + exp2 <- exp[, sd] + if (sum(is.na(exp2)) >= 1) { + app <- rep(NA, length(exp2)) + nas_pos <- which(is.na(exp2)) + exp2 <- exp2[!is.na(exp2)] + app[-nas_pos] <- doQmap(as.vector(exp2), adjust, ...) + } else { + app <- doQmap(as.vector(exp2), adjust, ...) + } + applied[, sd] <- app + } + } + } + } else { + applied <- exp_cor * NA + if (na.rm) { + tryCatch({ + adjust <- fitQmap(obs[!is.na(obs)], exp[!is.na(exp)], + method = method, ...) + applied[!is.na(exp_cor)] <- doQmap(exp_cor[!is.na(exp_cor)], + adjust, ...) + }, + error = function(error_message) { + return(applied) + }) + } else { + adjust <- fitQmap(as.vector(obs), as.vector(exp), method = method, ...) + applied <- doQmap(as.vector(exp_cor), adjust, ...) + } + dim(applied) <- dim(exp_cor) + } + return(applied) +}