diff --git a/DESCRIPTION b/DESCRIPTION index 7fe8d6812209852a7de75a7eb923886ee2847137..52999f2a9e329e9db763e3e2f5e482329bd4dc7b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,15 +1,17 @@ Package: CSTools Title: Assessing Skill of Climate Forecasts on Seasonal-to-Decadal Timescales -Version: 4.0.0 +Version: 4.0.1 Authors@R: c( person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-8568-3071")), person("Louis-Philippe", "Caron", , "louis-philippe.caron@bsc.es", role = "aut", comment = c(ORCID = "0000-0001-5221-0147")), person("Carmen", "Alvarez-Castro", , "carmen.alvarez-castro@cmcc.it", role = "aut", comment = c(ORCID = "0000-0002-9958-010X")), person("Lauriane", "Batte", , "lauriane.batte@meteo.fr", role = "aut"), + person("Carlos", "Delgado", , "carlos.delgado@bsc.es", role = "aut"), person("Jost", "von Hardenberg", , email = c("j.vonhardenberg@isac.cnr.it", "jost.hardenberg@polito.it"), role = "aut", comment = c(ORCID = "0000-0002-5312-8070")), person("Llorenç", "LLedo", , "llledo@bsc.es", role = "aut"), person("Nicolau", "Manubens", , "nicolau.manubens@bsc.es", role = "aut"), + person("Lluís", "Palma", , "lluis.palma@bsc.es", role = "aut"), person("Eroteida", "Sanchez-Garcia", , "esanchezg@aemet.es", role = "aut"), person("Bert", "van Schaeybroeck", , "bertvs@meteo.be", role = "aut"), person("Veronica", "Torralba", , "veronica.torralba@bsc.es", role = "aut"), diff --git a/NAMESPACE b/NAMESPACE index 7e313bc5e0cdbdfbd3f2735c6958745477e20276..34851f66cb1eef519700cde1b5194b2c3469ae31 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,6 +4,7 @@ export(AdamontAnalog) export(Analogs) export(BEI_PDFBest) export(BEI_Weights) +export(BiasCorrection) export(CST_Analogs) export(CST_AnalogsPredictors) export(CST_Anomaly) @@ -11,12 +12,14 @@ export(CST_BEI_Weighting) export(CST_BiasCorrection) export(CST_Calibration) export(CST_CategoricalEnsCombination) +export(CST_DynBiasCorrection) export(CST_EnsClustering) export(CST_Load) export(CST_MergeDims) export(CST_MultiEOF) export(CST_MultiMetric) export(CST_MultivarRMSE) +export(CST_ProxiesAttractor) export(CST_QuantileMapping) export(CST_RFSlope) export(CST_RFTemp) @@ -28,6 +31,7 @@ export(CST_SplitDim) export(CST_WeatherRegimes) export(Calibration) export(CategoricalEnsCombination) +export(DynBiasCorrection) export(EnsClustering) export(MergeDims) export(MultiEOF) @@ -37,6 +41,8 @@ export(PlotForecastPDF) export(PlotMostLikelyQuantileMap) export(PlotPDFsOLE) export(PlotTriangles4Categories) +export(Predictability) +export(ProxiesAttractor) export(QuantileMapping) export(RFSlope) export(RFTemp) @@ -114,4 +120,4 @@ importFrom(utils,head) importFrom(utils,read.table) importFrom(utils,tail) importFrom(verification,verify) -useDynLib(CSTools) +useDynLib(CSTools, .registration = TRUE) diff --git a/NEWS.md b/NEWS.md index 23d198894d3d5db5510465662479fd0a6947b591..6d66f0a29c3e7b22226d91d6ca505b8e9a462961 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,22 @@ +### CSTools 4.0.1 +**Submission date to CRAN: XX-06-2021** + +- New features: + + Dynamical Bias Correction method: `CST_ProxiesAttractors` and `CST_DynBiasCorrection` + (optionally `Predictability`) + + CST_BiasCorrection and BiasCorrection allows to calibrate a forecast given the calibration in the hindcast by using parameter 'exp_cor'. + + Use cases + + CST_SaveExp includes parameter extra_string + + PlotCombinedMap includes parameter cex_bar_titles + +- Fixes: + + Calibration retains correlation absolute value + + Calibration fixed when cal.methodi == rpc-based, apply_to == sign, + eval.method == 'leave-one-out' and the correlation is not significant + + PlotMostLikelyQuantileMap reoder latitudes of an array provided in 'dots' parameter. + ### CSTools 4.0.0 -**Submission date to CRAN: XX-12-2020** +**Submission date to CRAN: 23-02-2021** - New features: + ADAMONT downscaling method: requires CST_AdamontAnalogs and CST_AdamontQQCor functions diff --git a/R/CST_Analogs.R b/R/CST_Analogs.R index 24475de158c0d71f629d68e7b2448a50b5f9f802..1b463f4b9805f2c1606ed43cc991c04338c35bf9 100644 --- a/R/CST_Analogs.R +++ b/R/CST_Analogs.R @@ -170,7 +170,8 @@ CST_Analogs <- function(expL, obsL, expVar = NULL, obsVar = NULL, region = NULL, time_obsL <- obsL$Dates$start } res <- Analogs(expL$data, obsL$data, time_obsL = time_obsL, - time_expL = time_expL, expVar = expVar$data, + time_expL = time_expL, lonL = expL$lon, + latL = expL$lat, expVar = expVar$data, obsVar = obsVar$data, criteria = criteria, excludeTime = excludeTime, region = region, lonVar = as.vector(obsVar$lon), latVar = as.vector(obsVar$lat), @@ -253,6 +254,8 @@ CST_Analogs <- function(expL, obsL, expVar = NULL, obsVar = NULL, region = NULL, #'@param time_expL an array of N named dimensions (coinciding with time #'dimensions in expL) of character string(s) indicating the date(s) of the #'experiment in the format "dd/mm/yyyy". Time(s) to find the analogs. +#'@param lonL a vector containing the longitude of parameter 'expL'. +#'@param latL a vector containing the latitude of parameter 'expL'. #'@param excludeTime an array of N named dimensions (coinciding with time #'dimensions in expL) of character string(s) indicating the date(s) of the #'observations in the format "dd/mm/yyyy" to be excluded during the search of @@ -343,21 +346,21 @@ CST_Analogs <- function(expL, obsL, expVar = NULL, obsVar = NULL, region = NULL, #'region=c(lonmin = -1 ,lonmax = 2, latmin = 30, latmax = 33) #'Local_scale <- Analogs(expL = expSLP, obsL = obsSLP, time_obsL = time_obsSLP, #' obsVar = obs.pr, criteria = "Local_dist", -#' lonVar = seq(-1, 5, 1.5),latVar = seq(30, 35, 1.5), +#' lonL = seq(-1, 5, 1.5),latL = seq(30, 35, 1.5), #' region = region,time_expL = "01-10-2000", #' nAnalogs = 10, AnalogsInfo = TRUE) #' #'# Example 6: list of best analogs using criteria 'Local_dist' and 2 #'Local_scale <- Analogs(expL = expSLP, obsL = obsSLP, time_obsL = time_obsSLP, -#' criteria = "Local_dist", lonVar = seq(-1, 5, 1.5), -#' latVar = seq(30, 35, 1.5), region = region, +#' criteria = "Local_dist", lonL = seq(-1, 5, 1.5), +#' latL = seq(30, 35, 1.5), region = region, #' time_expL = "01-10-2000", nAnalogs = 5, #' AnalogsInfo = TRUE) #' #'# Example 7: Downscaling using Local_dist criteria #'Local_scale <- Analogs(expL = expSLP, obsL = obsSLP, time_obsL = time_obsSLP, -#' criteria = "Local_dist", lonVar = seq(-1, 5, 1.5), -#' latVar = seq(30, 35, 1.5), region = region, +#' criteria = "Local_dist", lonL = seq(-1, 5, 1.5), +#' latL = seq(30, 35, 1.5), region = region, #' time_expL = "01-10-2000", #' nAnalogs = 10, AnalogsInfo = FALSE) #' @@ -366,14 +369,16 @@ CST_Analogs <- function(expL, obsL, expVar = NULL, obsVar = NULL, region = NULL, #'dim(exp.pr) <- dim(expSLP) #'Local_scalecor <- Analogs(expL = expSLP, obsL = obsSLP,time_obsL = time_obsSLP, #' obsVar = obs.pr, expVar = exp.pr, -#' criteria = "Local_cor", lonVar = seq(-1, 5, 1.5), -#' time_expL = "01-10-2000", latVar = seq(30, 35, 1.5), +#' criteria = "Local_cor", lonL = seq(-1, 5, 1.5), +#' time_expL = "01-10-2000", latL = seq(30, 35, 1.5), +#' lonVar = seq(-1, 5, 1.5), latVar = seq(30, 35, 1.5), #' nAnalogs = 8, region = region, AnalogsInfo = FALSE) #'# same but without imposing nAnalogs,so nAnalogs will be set by default as 10 #'Local_scalecor <- Analogs(expL = expSLP, obsL = obsSLP,time_obsL = time_obsSLP, #' obsVar = obs.pr, expVar = exp.pr, -#' criteria = "Local_cor", lonVar = seq(-1,5,1.5), -#' time_expL = "01-10-2000", latVar=seq(30, 35, 1.5), +#' lonVar = seq(-1, 5, 1.5), latVar = seq(30, 35, 1.5), +#' criteria = "Local_cor", lonL = seq(-1,5,1.5), +#' time_expL = "01-10-2000", latL =seq(30, 35, 1.5), #' region = region, AnalogsInfo = TRUE) #' #'#'Example 9: List of best analogs in the three criterias Large_dist, @@ -382,11 +387,12 @@ CST_Analogs <- function(expL, obsL, expVar = NULL, obsVar = NULL, region = NULL, #' nAnalogs = 7, AnalogsInfo = TRUE) #'Local_scale <- Analogs(expL = expSLP, obsL = obsSLP, time_obsL = time_obsSLP, #' time_expL = "01-10-2000", criteria = "Local_dist", -#' lonVar = seq(-1, 5, 1.5), latVar = seq(30, 35, 1.5), +#' lonL = seq(-1, 5, 1.5), latL = seq(30, 35, 1.5), #' nAnalogs = 7,region = region, AnalogsInfo = TRUE) #'Local_scalecor <- Analogs(expL = expSLP, obsL = obsSLP,time_obsL = time_obsSLP, #' obsVar = obsSLP, expVar = expSLP, #' time_expL = "01-10-2000",criteria = "Local_cor", +#' lonL = seq(-1, 5, 1.5), latL = seq(30, 35, 1.5), #' lonVar = seq(-1, 5, 1.5), latVar = seq(30, 35, 1.5), #' nAnalogs = 7,region = region, #' AnalogsInfo = TRUE) @@ -404,7 +410,8 @@ CST_Analogs <- function(expL, obsL, expVar = NULL, obsVar = NULL, region = NULL, #' time_obsL = time_obsSLP, time_expL = time_expSLP, #' excludeTime = excludeTime, AnalogsInfo = TRUE) #'@export -Analogs <- function(expL, obsL, time_obsL, time_expL = NULL, expVar = NULL, +Analogs <- function(expL, obsL, time_obsL, time_expL = NULL, + lonL = NULL, latL = NULL, expVar = NULL, obsVar = NULL, criteria = "Large_dist",excludeTime = NULL, lonVar = NULL, latVar = NULL, region = NULL, @@ -640,6 +647,7 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL, expVar = NULL, fun = .analogs, time_obsL, expVar = expVar, time_expL=time_expL, excludeTime=excludeTime, obsVar = obsVar, criteria = criteria, + lonL = lonL, latL = latL, lonVar = lonVar, latVar = latVar, region = region, nAnalogs = nAnalogs, AnalogsInfo = AnalogsInfo, output_dims = c('nAnalogs', 'lat', 'lon'), @@ -652,6 +660,7 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL, expVar = NULL, fun = .analogs,time_obsL, time_expL=time_expL, excludeTime=excludeTime, expVar = expVar, criteria = criteria, + lonL = lonL, latL = latL, lonVar = lonVar, latVar = latVar, region = region, nAnalogs = nAnalogs, AnalogsInfo = AnalogsInfo, output_dims = c('nAnalogs', 'lat', 'lon'), @@ -664,6 +673,7 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL, expVar = NULL, fun = .analogs, criteria = criteria,time_obsL, time_expL=time_expL, excludeTime=excludeTime, + lonL = lonL, latL = latL, lonVar = lonVar, latVar = latVar, region = region, nAnalogs = nAnalogs, AnalogsInfo = AnalogsInfo, output_dims = c('nAnalogs', 'lat', 'lon'), @@ -676,6 +686,7 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL, expVar = NULL, fun = .analogs, time_obsL, expVar = expVar, time_expL=time_expL, excludeTime=excludeTime, obsVar = obsVar, criteria = criteria, + lonL = lonL, latL = latL, lonVar = lonVar, latVar = latVar, region = region, nAnalogs = nAnalogs, AnalogsInfo = AnalogsInfo, output_dims = list(fields = c('nAnalogs', 'lat', 'lon'), @@ -690,6 +701,7 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL, expVar = NULL, fun = .analogs,time_obsL, time_expL=time_expL, excludeTime=excludeTime, expVar = expVar, criteria = criteria, + lonL = lonL, latL = latL, lonVar = lonVar, latVar = latVar, region = region, nAnalogs = nAnalogs, AnalogsInfo = AnalogsInfo, output_dims = list(fields = c('nAnalogs', 'lat', 'lon'), @@ -705,6 +717,7 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL, expVar = NULL, fun = .analogs,time_obsL, criteria = criteria, time_expL=time_expL, excludeTime=excludeTime, + lonL = lonL, latL = latL, lonVar = lonVar, latVar = latVar, region = region, nAnalogs = nAnalogs, AnalogsInfo = AnalogsInfo, output_dims = list(fields = c('nAnalogs', 'lat', 'lon'), @@ -719,6 +732,7 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL, expVar = NULL, .analogs <- function(expL, obsL, time_expL, excludeTime = NULL, obsVar = NULL, expVar = NULL, time_obsL, criteria = "Large_dist", + lonL = NULL, latL = NULL, lonVar = NULL, latVar = NULL, region = NULL, nAnalogs = NULL, AnalogsInfo = FALSE) { @@ -796,7 +810,8 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL, expVar = NULL, expVar = expVar, obsVar = obsVar, criteria = criteria, AnalogsInfo = AnalogsInfo, - nAnalogs = nAnalogs,lonVar = lonVar, + nAnalogs = nAnalogs, + lonL = lonL, latL = latL, lonVar = lonVar, latVar = latVar, region = region) if (AnalogsInfo == TRUE) { return(list(AnalogsFields = Analog_result$AnalogsFields, @@ -807,14 +822,17 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL, expVar = NULL, return(AnalogsFields = Analog_result$AnalogsFields) } } -FindAnalog <- function(expL, obsL, time_obsL, expVar, obsVar, criteria, lonVar, +FindAnalog <- function(expL, obsL, time_obsL, expVar, obsVar, criteria, + lonL, latL, lonVar, latVar, region, nAnalogs = nAnalogs, AnalogsInfo = AnalogsInfo) { position <- Select(expL = expL, obsL = obsL, expVar = expVar, - obsVar = obsVar, criteria = criteria, lonVar = lonVar, + obsVar = obsVar, criteria = criteria, + lonL = lonL, latL = latL, lonVar = lonVar, latVar = latVar, region = region)$position metrics<- Select(expL = expL, obsL = obsL, expVar = expVar, - obsVar = obsVar, criteria = criteria, lonVar = lonVar, + obsVar = obsVar, criteria = criteria, lonL = lonL, + latL = latL, lonVar = lonVar, latVar = latVar, region = region)$metric.original best <- Apply(list(position), target_dims = c('time', 'pos'), fun = BestAnalog, criteria = criteria, @@ -824,8 +842,8 @@ FindAnalog <- function(expL, obsL, time_obsL, expVar, obsVar, criteria, lonVar, dim(Analogs_dates) <- dim(best) if (all(!is.null(region), !is.null(lonVar), !is.null(latVar))) { if (is.null(obsVar)) { - obsVar <- SelBox(obsL, lon = lonVar, lat = latVar, region = region)$data - expVar <- SelBox(expL, lon = lonVar, lat = latVar, region=region)$data + obsVar <- SelBox(obsL, lon = lonL, lat = latL, region = region)$data + expVar <- SelBox(expL, lon = lonL, lat = latL, region=region)$data Analogs_fields <- Subset(obsVar, along = which(names(dim(obsVar)) == 'time'), indices = best) @@ -955,7 +973,7 @@ BestAnalog <- function(position, nAnalogs = nAnalogs, AnalogsInfo = FALSE, } } Select <- function(expL, obsL, expVar = NULL, obsVar = NULL, - criteria = "Large_dist", + criteria = "Large_dist", lonL = NULL, latL = NULL, lonVar = NULL, latVar = NULL, region = NULL) { names(dim(expL)) <- replace_repeat_dimnames(names(dim(expL)), names(dim(obsL))) @@ -989,11 +1007,11 @@ Select <- function(expL, obsL, expVar = NULL, obsVar = NULL, position = pos1)) } if (criteria == "Local_dist" | criteria == "Local_cor") { - obs <- SelBox(obsL, lon = lonVar, lat = latVar, region = region)$data - exp <- SelBox(expL, lon = lonVar, lat = latVar, region = region)$data + obs <- SelBox(obsL, lon = lonL, lat = latL, region = region)$data + exp <- SelBox(expL, lon = lonL, lat = latL, region = region)$data metric2 <- Apply(list(obs), target_dims = list(c('lat', 'lon')), fun = .select, exp, metric = "dist")$output1 - metric2.original=metric2 + metric2.original = metric2 dim(metric2) <- c(dim(metric2), metric=1) margins <- c(1 : (length(dim(metric2))))[-dim_time_obs] pos2 <- apply(metric2, margins, order) diff --git a/R/CST_BiasCorrection.R b/R/CST_BiasCorrection.R index 9baf897bb680eff59d5b3bd4b3169e4a880ee65b..0b8375130cb170345fd7e356214f577b3703662b 100644 --- a/R/CST_BiasCorrection.R +++ b/R/CST_BiasCorrection.R @@ -5,6 +5,7 @@ #' #'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data} #'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed data in the element named \code{$data}. +#'@param exp_cor an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonl forecast experiment to be corrected. If it is NULL, the 'exp' forecast will be corrected. #'@param na.rm a logical value indicating whether missing values should be stripped before the computation proceeds, by default it is set to FALSE. #' #'@return an object of class \code{s2dv_cube} containing the bias corrected forecasts in the element called \code{$data} with the same dimensions of the experimental data. @@ -30,28 +31,70 @@ #'a <- CST_BiasCorrection(exp = exp, obs = obs) #'str(a) #'@export -CST_BiasCorrection <- function(exp, obs, na.rm = FALSE) { +CST_BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE) { 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 (dim(obs$data)['member'] != 1) { stop("The length of the dimension 'member' in the component 'data' ", "of the parameter 'obs' must be equal to 1.") } - dimnames <- names(dim(exp$data)) - BiasCorrected <- BiasCorrection(exp = exp$data, obs = obs$data, na.rm = na.rm) - pos <- match(dimnames, names(dim(BiasCorrected))) - BiasCorrected <- aperm(BiasCorrected, pos) - names(dim(BiasCorrected)) <- dimnames - exp$data <- BiasCorrected - exp$Datasets <- c(exp$Datasets, obs$Datasets) - exp$source_files <- c(exp$source_files, obs$source_files) - return(exp) + if (is.null(exp_cor)) { + dimnames <- names(dim(exp$data)) + BiasCorrected <- BiasCorrection(exp = exp$data, obs = obs$data, na.rm = na.rm) + pos <- match(dimnames, names(dim(BiasCorrected))) + BiasCorrected <- aperm(BiasCorrected, pos) + names(dim(BiasCorrected)) <- dimnames + exp$data <- BiasCorrected + exp$Datasets <- c(exp$Datasets, obs$Datasets) + exp$source_files <- c(exp$source_files, obs$source_files) + return(exp) + } else { + if (!inherits(exp_cor, 's2dv_cube')) { + stop("Parameter 'var_cor' must be of the class 's2dv_cube'.") + } + dimnames <- names(dim(exp_cor$data)) + BiasCorrected <- BiasCorrection(exp = exp$data, obs = obs$data, + exp_cor = exp_cor$data, na.rm = na.rm) + pos <- match(dimnames, names(dim(BiasCorrected))) + BiasCorrected <- aperm(BiasCorrected, pos) + names(dim(BiasCorrected)) <- dimnames + exp_cor$data <- BiasCorrected + 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) + } } - -BiasCorrection <- function (exp, obs , na.rm = FALSE) { +#' Bias Correction based on the mean and standard deviation adjustment +#' +#'@author Verónica Torralba, \email{veronica.torralba@bsc.es} +#'@description This function applies the simple bias adjustment technique described in Torralba et al. (2017). The adjusted forecasts have an equivalent standard deviation and mean to that of the reference dataset. +#' +#'@param exp a multidimensional array with named dimensions containing the seasonal forecast experiment data with at least 'member' and 'sdate' dimensions. +#'@param obs a multidimensional array with named dimensions containing the observed data with at least 'sdate' dimension. +#'@param exp_cor a multidimensional array with named dimensions containing the seasonl forecast experiment to be corrected. If it is NULL, the 'exp' forecast will be corrected. +#'@param na.rm a logical value indicating whether missing values should be stripped before the computation proceeds, by default it is set to FALSE. +#' +#'@return an object of class \code{s2dv_cube} containing the bias corrected forecasts in the element called \code{$data} with the same dimensions of the experimental data. +#' +#'@references Torralba, V., F.J. Doblas-Reyes, D. MacLeod, I. Christel and M. Davis (2017). Seasonal climate prediction: a new source of information for the management of wind energy resources. Journal of Applied Meteorology and Climatology, 56, 1231-1247, doi:10.1175/JAMC-D-16-0204.1. (CLIM4ENERGY, EUPORIAS, NEWA, RESILIENCE, SPECS) +#' +#'@import multiApply +#'@examples +#' +#'# Example +#'# Creation of sample s2dverification objects. These are not complete +#'# s2dverification objects though. The Load function returns complete objects. +#'mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7) +#'dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7) +#'obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7) +#'dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7) +#'a <- BiasCorrection(exp = mod1, obs = obs1) +#'str(a) +#'@export +BiasCorrection <- function (exp, obs, exp_cor = NULL, na.rm = FALSE) { if (!all(c('member', 'sdate') %in% names(dim(exp)))) { stop("Parameter 'exp' must have the dimensions 'member' and 'sdate'.") @@ -83,37 +126,55 @@ BiasCorrection <- function (exp, obs , na.rm = FALSE) { if ('member' %in% names(dim(obs))) { target_dims_obs <- c('member', target_dims_obs) } - - BiasCorrected <- Apply(data = list(var_obs = obs, var_exp = exp), - target_dims = list(target_dims_obs, c('member', 'sdate')), + if (is.null(exp_cor)) { + BiasCorrected <- Apply(data = list(var_obs = obs, var_exp = exp), + target_dims = list(target_dims_obs, + c('member', 'sdate')), fun = .sbc , na.rm = na.rm)$output1 + } else { + BiasCorrected <- Apply(data = list(var_obs = obs, + var_exp = exp, + var_cor = exp_cor), + target_dims = list(target_dims_obs, + c('member', 'sdate'), + c('member', 'sdate')), + fun = .sbc , output_dims = c('member','sdate'), + na.rm = na.rm)$output1 + } return(BiasCorrected) } -.sbc <- function(var_obs, var_exp , na.rm = FALSE) { +.sbc <- function(var_obs, var_exp , var_cor = NULL, na.rm = FALSE) { nmembers <- dim(var_exp)['member'][] ntime <- dim(var_exp)['sdate'][] - if (all(names(dim(var_exp)) != c('member','sdate'))) { - var_exp <- t(var_exp) - } + #if (all(names(dim(var_exp)) != c('member','sdate'))) { + # var_exp <- t(var_exp) + #} corrected <- NA * var_exp - - for (t in 1 : ntime) { - # defining forecast,hindcast and observation in cross-validation - fcst <- var_exp[ , t] - hcst <- var_exp[ , -t] - obs <- var_obs[-t] - + + if (is.null(var_cor)){ + for (t in 1 : ntime) { + # parameters + sd_obs <- sd(var_obs[-t], na.rm = na.rm) + sd_exp <- sd(var_exp[ , -t], na.rm = na.rm) + clim_exp <- mean(var_exp[ , -t], na.rm = na.rm) + clim_obs <- mean(var_obs[-t], na.rm = na.rm) + + # bias corrected forecast + corrected[ , t] <- ((var_exp[ , t] - clim_exp) * (sd_obs / sd_exp)) + clim_obs + names(dim(corrected)) <- c('member', 'sdate') + } + } else { # parameters - sd_obs <- sd(obs , na.rm = na.rm) - sd_exp <- sd(hcst , na.rm = na.rm) - clim_exp <- mean(hcst , na.rm = na.rm) - clim_obs <- mean(obs , na.rm = na.rm) + sd_obs <- sd(var_obs, na.rm = na.rm) + sd_exp <- sd(var_exp, na.rm = na.rm) + clim_exp <- mean(var_exp, na.rm = na.rm) + clim_obs <- mean(var_obs, na.rm = na.rm) # bias corrected forecast - corrected[ , t] <- ((fcst - clim_exp) * (sd_obs / sd_exp)) + clim_obs + corrected <- ((var_cor - clim_exp) * (sd_obs / sd_exp)) + clim_obs + names(dim(corrected)) <- c('member') } - names(dim(corrected)) <- c('member', 'sdate') return(corrected) } diff --git a/R/CST_Calibration.R b/R/CST_Calibration.R index 58ef720777bd5ebb58cb5a3916b48ccbe71a6414..2884712ac89a94ba3203b5102c92e3c61afa1e05 100644 --- a/R/CST_Calibration.R +++ b/R/CST_Calibration.R @@ -72,7 +72,7 @@ CST_Calibration <- function(exp, obs, cal.method = "mse_min", #' #'@author Verónica Torralba, \email{veronica.torralba@bsc.es} #'@author Bert Van Schaeybroeck, \email{bertvs@meteo.be} -#'@description Four types of member-by-member bias correction can be performed. The \code{"bias"} method corrects the bias only, the \code{"evmos"} method applies a variance inflation technique to ensure the correction of the bias and the correspondence of variance between forecast and observation (Van Schaeybroeck and Vannitsem, 2011). The ensemble calibration methods \code{"mse_min"} and \code{"crps_min"} correct the bias, the overall forecast variance and the ensemble spread as described in Doblas-Reyes et al. (2005) and Van Schaeybroeck and Vannitsem (2015), respectively. While the \code{"mse_min"} method minimizes a constrained mean-squared error using three parameters, the \code{"crps_min"} method features four parameters and minimizes the Continuous Ranked Probability Score (CRPS). The \code{"rpc-based"} method adjusts the forecast variance ensuring that the ratio of predictable components (RPC) is equal to one, as in Eade et al. (2014). +#'@description Five types of member-by-member bias correction can be performed. The \code{"bias"} method corrects the bias only, the \code{"evmos"} method applies a variance inflation technique to ensure the correction of the bias and the correspondence of variance between forecast and observation (Van Schaeybroeck and Vannitsem, 2011). The ensemble calibration methods \code{"mse_min"} and \code{"crps_min"} correct the bias, the overall forecast variance and the ensemble spread as described in Doblas-Reyes et al. (2005) and Van Schaeybroeck and Vannitsem (2015), respectively. While the \code{"mse_min"} method minimizes a constrained mean-squared error using three parameters, the \code{"crps_min"} method features four parameters and minimizes the Continuous Ranked Probability Score (CRPS). The \code{"rpc-based"} method adjusts the forecast variance ensuring that the ratio of predictable components (RPC) is equal to one, as in Eade et al. (2014). #'@description Both in-sample or our out-of-sample (leave-one-out cross validation) calibration are possible. #'@references Doblas-Reyes F.J, Hagedorn R, Palmer T.N. The rationale behind the success of multi-model ensembles in seasonal forecasting-II calibration and combination. Tellus A. 2005;57:234-252. doi:10.1111/j.1600-0870.2005.00104.x #'@references Eade, R., Smith, D., Scaife, A., Wallace, E., Dunstone, N., Hermanson, L., & Robinson, N. (2014). Do seasonal-to-decadal climate predictions underestimate the predictability of the read world? Geophysical Research Letters, 41(15), 5620-5628. doi: 10.1002/2014GL061146 @@ -300,8 +300,8 @@ Calibration <- function(exp, obs, cal.method = "mse_min", #correct evaluation subset var.cor.fc[ , eval.dexes] <- .correct.crps.min.fc(fc.ev , mbm.par, na.rm = na.rm) } else if (cal.method == 'rpc-based') { - ens_mean.ev <- s2dv::MeanDims(data = fc.ev, dims = names(amt.mbr), na.rm = na.rm) - ens_mean.tr <- s2dv::MeanDims(data = fc.tr, dims = names(amt.mbr), na.rm = na.rm) ## Ensemble mean + ens_mean.ev <- multiApply::Apply(data = fc.ev, target_dims = names(amt.mbr), fun = mean, na.rm = na.rm)$output1 ## Ensemble mean + ens_mean.tr <- multiApply::Apply(data = fc.tr, target_dims = names(amt.mbr), fun = mean, na.rm = na.rm)$output1 ## Ensemble mean ens_spread.tr <- multiApply::Apply(data = list(fc.tr, ens_mean.tr), target_dims = names(amt.sdate), fun = "-")$output1 ## Ensemble spread exp_mean.tr <- mean(fc.tr, na.rm = na.rm) ## Mean (climatology) var_signal.tr <- var(ens_mean.tr, na.rm = na.rm) ## Ensemble mean variance @@ -314,7 +314,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", order = names(dims.fc)) dim(var.cor.fc) <- dims.fc } else { ## no significant -> replacing with observed climatology - var.cor.fc[ , eval.dexes] <- array(data = mean(obs.tr, na.rm = na.rm), dim = dim(fc.tr)) + var.cor.fc[ , eval.dexes] <- array(data = mean(obs.tr, na.rm = na.rm), dim = dim(fc.ev)) } } else { stop("unknown calibration method: ",cal.method) @@ -416,7 +416,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", } else { par.out[3] <- with(quant.obs.fc, obs.sd * sqrt(1. - cor.obs.fc^2) / fc.dev.sd) } - par.out[2] <- with(quant.obs.fc, cor.obs.fc * obs.sd / fc.ens.av.sd) + par.out[2] <- with(quant.obs.fc, abs(cor.obs.fc) * obs.sd / fc.ens.av.sd) par.out[1] <- with(quant.obs.fc, obs.av - par.out[2] * fc.ens.av.av, na.rm = na.rm) return(par.out) @@ -473,4 +473,4 @@ Calibration <- function(exp, obs, cal.method = "mse_min", .CalibrationMembersRPC <- function(exp, ens_mean, ens_mean_cal, var_obs, var_noise, r){ member_cal <- (exp - ens_mean) * sqrt(var_obs) * sqrt(1 - r^2) / sqrt(var_noise) + ens_mean_cal return(member_cal) -} \ No newline at end of file +} diff --git a/R/CST_DynBiasCorrection.R b/R/CST_DynBiasCorrection.R new file mode 100644 index 0000000000000000000000000000000000000000..20c263c6732bb0148a7514b482199f95d05af01a --- /dev/null +++ b/R/CST_DynBiasCorrection.R @@ -0,0 +1,272 @@ +#'@rdname CST_DynBiasCorrection +#'@title Performing a Bias Correction conditioned by the dynamical +#'properties of the data. +#' +#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} +#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} +#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it} +#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +#' +#'@description This function perform a bias correction conditioned by the +#'dynamical properties of the dataset. This function internally uses the functions +#''Predictability' to divide in terciles the two dynamical proxies +#'computed with 'CST_ProxiesAttractor'. A bias correction +#'between the model and the observations is performed using the division into +#'terciles of the local dimension 'dim' and inverse of the persistence 'theta'. +#'For instance, model values with lower 'dim' will be corrected with observed +#'values with lower 'dim', and the same for theta. The function gives two options +#'of bias correction: one for 'dim' and/or one for 'theta' +#' +#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., +#'and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large +#'scale atmospheric predictability.Nature Communications, 10(1), 1316. +#'DOI = https://doi.org/10.1038/s41467-019-09305-8 " +#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +#' Dynamical proxies of North Atlantic predictability and extremes. +#' Scientific Reports, 7-41278, 2017. +#' +#'@param exp an s2v_cube object with the experiment data +#'@param obs an s2dv_cube object with the reference data +#'@param method a character string indicating the method to apply bias +#'correction among these ones: "PTF","RQUANT","QUANT","SSPLIN" +#'@param wetday logical indicating whether to perform wet day correction +#'or not OR a numeric threshold below which all values are set to zero (by +#'default is set to 'FALSE'). +#'@param proxy a character string indicating the proxy for local dimension +#' 'dim' or inverse of persistence 'theta' to apply the dynamical +#' conditioned bias correction method. +#'@param quanti a number lower than 1 indicating the quantile to perform +#'the computation of local dimension and theta +#'@param ncores The number of cores to use in parallel computation +#' +#'@return dynbias an s2dvcube object with a bias correction performed +#'conditioned by local dimension 'dim' or inverse of persistence 'theta' +#' +#'@examples +#'# example 1: simple data s2dvcube style +#' set.seed(1) +#' expL <- rnorm(1:2000) +#' dim (expL) <- c(time =100,lat = 4, lon = 5) +#' obsL <- c(rnorm(1:1980),expL[1,,]*1.2) +#' dim (obsL) <- c(time = 100,lat = 4, lon = 5) +#' time_obsL <- paste(rep("01", 100), rep("01", 100), 1920 : 2019, sep = "-") +#' time_expL <- paste(rep("01", 100), rep("01", 100), 1929 : 2019, sep = "-") +#' lon <- seq(-1,5,1.5) +#' lat <- seq(30,35,1.5) +#' # qm=0.98 # too high for this short dataset, it is possible that doesn't +#' # get the requirement, in that case it would be necessary select a lower qm +#' # for instance qm=0.60 +#' expL <- s2dv_cube(data = expL, lat = lat, lon = lon, +#' Dates = list(start = time_expL, end = time_expL)) +#' obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, +#' Dates = list(start = time_obsL, end = time_obsL)) +#' # to use DynBiasCorrection +#' dynbias1 <- DynBiasCorrection(exp = expL$data, obs = obsL$data, proxy= "dim", +#' quanti = 0.6) +#' # to use CST_DynBiasCorrection +#' dynbias2 <- CST_DynBiasCorrection(exp = expL, obs = obsL, proxy= "dim", +#' quanti = 0.6) +#' +#'@export +CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', wetday=FALSE, + proxy = "dim", quanti, + ncores = NULL) { + if (!inherits(obs, 's2dv_cube')) { + stop("Parameter 'obs' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + if (!inherits(exp, 's2dv_cube')) { + stop("Parameter 'exp' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + exp$data <- DynBiasCorrection(exp = exp$data, obs = obs$data, method = method, + wetday=wetday, + proxy = proxy, quanti = quanti, ncores = ncores) + return(exp) +} +#'@rdname DynBiasCorrection +#'@title Performing a Bias Correction conditioned by the dynamical +#'properties of the data. +#' +#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} +#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} +#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it} +#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +#' +#'@description This function perform a bias correction conditioned by the +#'dynamical properties of the dataset. This function used the functions +#''CST_Predictability' to divide in terciles the two dynamical proxies +#'computed with 'CST_ProxiesAttractor'. A bias correction +#'between the model and the observations is performed using the division into +#'terciles of the local dimension 'dim' and inverse of the persistence 'theta'. +#'For instance, model values with lower 'dim' will be corrected with observed +#'values with lower 'dim', and the same for theta. The function gives two options +#'of bias correction: one for 'dim' and/or one for 'theta' +#' +#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., +#'and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large +#'scale atmospheric predictability.Nature Communications, 10(1), 1316. +#'DOI = https://doi.org/10.1038/s41467-019-09305-8 " +#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +#' Dynamical proxies of North Atlantic predictability and extremes. +#' Scientific Reports, 7-41278, 2017. +#' +#'@param exp a multidimensional array with named dimensions with the +#'experiment data +#'@param obs a multidimensional array with named dimensions with the +#'observation data +#'@param method a character string indicating the method to apply bias +#'correction among these ones: +#'"PTF","RQUANT","QUANT","SSPLIN" +#'@param wetday logical indicating whether to perform wet day correction +#'or not OR a numeric threshold below which all values are set to zero (by +#'default is set to 'FALSE'). +#'@param proxy a character string indicating the proxy for local dimension +#''dim' or inverse of persistence 'theta' to apply the dynamical conditioned +#'bias correction method. +#'@param quanti a number lower than 1 indicating the quantile to perform the +#'computation of local dimension and theta +#'@param ncores The number of cores to use in parallel computation +#' +#'@return a multidimensional array with named dimensions with a bias correction +#'performed conditioned by local dimension 'dim' or inverse of persistence 'theta' +#' +#'@import multiApply +#'@importFrom s2dverification Subset +#'@import qmap +#'@examples +#'expL <- rnorm(1:2000) +#'dim (expL) <- c(time =100,lat = 4, lon = 5) +#'obsL <- c(rnorm(1:1980),expL[1,,]*1.2) +#'dim (obsL) <- c(time = 100,lat = 4, lon = 5) +#'dynbias <- DynBiasCorrection(exp = expL, obs = obsL, method='QUANT', +#' proxy= "dim", quanti = 0.6) +#'@export +DynBiasCorrection<- function(exp, obs, method = 'QUANT',wetday=FALSE, + proxy = "dim", quanti, ncores = NULL){ + if (is.null(obs)) { + stop("Parameter 'obs' cannot be NULL.") + } + if (is.null(exp)) { + stop("Parameter 'exp' cannot be NULL.") + } + if (is.null(method)) { + stop("Parameter 'method' cannot be NULL.") + } + if (is.null(quanti)) { + stop("Parameter 'quanti' cannot be NULL.") + } + if (is.null(proxy)) { + stop("Parameter 'proxy' cannot be NULL.") + } + dims <- dim(exp) + + attractor.obs <- ProxiesAttractor(data = obs, quanti = quanti) + predyn.obs <- Predictability(dim = attractor.obs$dim, + theta = attractor.obs$theta) + attractor.exp <- ProxiesAttractor(data = exp, quanti = quanti) + predyn.exp <- Predictability(dim = attractor.exp$dim, + theta = attractor.exp$theta) + + if (!(any(names(dim(exp)) %in% 'time'))){ + if (any(names(dim(exp)) %in% 'sdate')) { + if (any(names(dim(exp)) %in% 'ftime')) { + exp <- MergeDims(exp, merge_dims = c('ftime', 'sdate'), + rename_dim = 'time') + } + } + } + if (!(any(names(dim(obs)) %in% 'time'))){ + if (any(names(dim(obs)) %in% 'sdate')) { + if (any(names(dim(obs)) %in% 'ftime')) { + obs <- MergeDims(obs, merge_dims = c('ftime', 'sdate'), + rename_dim = 'time') + } + } + } + + dim_exp <- dim(exp) + names_to_check <- names(dim_exp)[which(names(dim_exp) %in% + c('time', 'lat', 'lon', 'sdate') == FALSE)] + if (length(names_to_check) > 0) { + dim_obs <- dim(obs) + if (any(names(dim_obs) %in% names_to_check)) { + if (any(dim_obs[which(names(dim_obs) %in% names_to_check)] != + dim_exp[which(names(dim_exp) %in% names_to_check)])) { + for (i in names_to_check) { + pos <- which(names(dim_obs) == i) + names(dim(obs))[pos] <- ifelse(dim_obs[pos] != + dim_exp[which(names(dim_exp) == i)], + paste0('obs_', names(dim_obs[pos])), + names(dim(obs)[pos])) + } + warning("Common dimension names with different length are renamed.") + } + } + } + + if (proxy == "dim") { + adjusted <- Apply(list(exp, obs), target_dims = 'time', + fun = .dynbias, method, wetday, + predyn.exp = predyn.exp$pred.dim$pos.d, + predyn.obs = predyn.obs$pred.dim$pos.d, + ncores = ncores, output_dims = 'time')$output1 + } else if (proxy == "theta") { + adjusted <- Apply(list(exp, obs), target_dims = 'time', + fun = .dynbias, method, wetday, + predyn.exp = predyn.exp$pred.theta$pos.t, + predyn.obs = predyn.obs$pred.theta$pos.t, + ncores = ncores, output_dims = 'time')$output1 + } else { + stop ("Parameter 'proxy' must be set as 'dim' or 'theta'.") + } + + if(any(names(dim(adjusted)) %in% 'memberObs')){ + if(dim(adjusted)['memberObs'] == 1){ + adjusted <- Subset(adjusted,along='memberObs',indices=1,drop = 'selected') + }else{ + print('Dimension member in obs changed to memberObs') + } + } + + if(any(names(dim(adjusted)) %in% 'datasetObs')){ + if(dim(adjusted)['datasetObs'] == 1){ + adjusted <- Subset(adjusted,along='datasetObs',indices=1,drop = 'selected') + }else{ + print('Dimension dataset in obs changed to datasetObs') + } + } + return(adjusted) +} + +.dynbias <- function(exp, obs, method, wetday,predyn.exp, predyn.obs) { + result <- array(rep(NA, length(exp))) + res <- lapply(1:3, function(x) { + exp_sub <- exp[predyn.exp[[x]]] + obs_sub <- obs[predyn.obs[[x]]] + adjust <- .qbiascorrection(exp_sub, obs_sub, method,wetday) + result[predyn.exp[[x]]] <<- adjust + return(NULL) + }) + return(result) +} +.qbiascorrection <- function(expX, obsX, method,wetday) { + ## functions fitQmap and doQmap + if (method == "PTF") { + qm.fit <- fitQmap(obsX, expX, method = "PTF", transfun = "expasympt", + cost = "RSS", wet.day = wetday) + qmap <- doQmap(expX, qm.fit) + } else if (method == "QUANT") { + qm.fit <- fitQmap(obsX, expX, method = "QUANT", qstep = 0.01, wet.day = wetday) + qmap <- doQmap(expX, qm.fit, type = "tricub") + } else if (method == "RQUANT") { + qm.fit <- fitQmap(obsX, expX, method = "RQUANT", qstep = 0.01,wet.day = wetday) + qmap <- doQmap(expX, qm.fit, type = "linear") + } else if (method == "SSPLIN") { + qm.fit <- fitQmap(obsX, expX, qstep = 0.01, method = "SSPLIN",wet.day = wetday) + qmap <- doQmap(expX, qm.fit) + } else { + stop ("Parameter 'method' doesn't match any of the available methods.") + } + return(qmap) +} diff --git a/R/CST_ProxiesAttractor.R b/R/CST_ProxiesAttractor.R new file mode 100644 index 0000000000000000000000000000000000000000..518968fa0bbd53b912beb01efee4ccc233779d06 --- /dev/null +++ b/R/CST_ProxiesAttractor.R @@ -0,0 +1,180 @@ +#'@rdname CST_ProxiesAttractor +#'@title Computing two dinamical proxies of the attractor in s2dv_cube. +#' +#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} +#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} +#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it} +#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +#' +#'@description This function computes two dinamical proxies of the attractor: +#'The local dimension (d) and the inverse of the persistence (theta) for an +#''s2dv_cube' object. +#'These two parameters will be used as a condition for the computation of dynamical +#'scores to measure predictability and to compute bias correction conditioned by +#'the dynamics with the function DynBiasCorrection +#'Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in +#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019). +#' The hammam effect or how a warm ocean enhances large scale atmospheric predictability. +#' Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 " +#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +#' Dynamical proxies of North Atlantic predictability and extremes. +#' Scientific Reports, 7-41278, 2017. +#' +#'@param data a s2dv_cube object with the data to create the attractor. Must be a matrix with the timesteps in nrow +#'and the grids in ncol(dat(time,grids) +# +#'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta +#' +#'@param ncores The number of cores to use in parallel computation +#' +#'@return dim and theta +#' +#'@examples +#'# Example 1: Computing the attractor using simple s2dv data +#'attractor <- CST_ProxiesAttractor(data = lonlat_data$obs, quanti = 0.6) +#' +#'@export +CST_ProxiesAttractor <- function(data, quanti, ncores = NULL){ + if (!inherits(data, 's2dv_cube')) { + stop("Parameter 'data' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + if (is.null(quanti)) { + stop("Parameter 'quanti' cannot be NULL.") + } + + data$data <- ProxiesAttractor(data = data$data, quanti = quanti, ncores = ncores) + + return(data) +} +#'@rdname ProxiesAttractor +#' +#'@title Computing two dinamical proxies of the attractor. +#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} +#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} +#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it} +#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +#' +#'@description This function computes two dinamical proxies of the attractor: +#'The local dimension (d) and the inverse of the persistence (theta). +#'These two parameters will be used as a condition for the computation of dynamical +#'scores to measure predictability and to compute bias correction conditioned by +#'the dynamics with the function DynBiasCorrection. +#'Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in: +#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019). +#' The hammam effect or how a warm ocean enhances large scale atmospheric predictability. +#' Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 " +#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +#' Dynamical proxies of North Atlantic predictability and extremes. +#' Scientific Reports, 7-41278, 2017. +#' +#'@param data a multidimensional array with named dimensions to create the attractor. It requires a temporal dimension named 'time' and spatial dimensions called 'lat' and 'lon', or 'latitude' and 'longitude' or 'grid'. +#'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta +#'@param ncores The number of cores to use in parallel computation. +#' +#'@return dim and theta +#' +#'@import multiApply +#' +#'@examples +#'# Example 1: Computing the attractor using simple data +#'# Creating an example of matrix data(time,grids): +#'mat <- array(rnorm(36 * 40), c(time = 36, grid = 40)) +#'qm <- 0.90 # imposing a threshold +#'Attractor <- ProxiesAttractor(data = mat, quanti = qm) +#'# to plot the result +#'time = c(1:length(Attractor$theta)) +#'layout(matrix(c(1, 3, 2, 3), 2, 2)) +#'plot(time, Attractor$dim, xlab = 'time', ylab = 'd', +#' main = 'local dimension', type = 'l') +#'plot(time, Attractor$theta, xlab = 'time', ylab = 'theta', main = 'theta') +#'plot(Attractor$dim, Attractor$theta, col = 'blue', +#' main = "Proxies of the Attractor", +#' xlab = "local dimension", ylab = "theta", lwd = 8, 'p') +#' +#'@export + +ProxiesAttractor <- function(data, quanti, ncores = NULL){ + if (is.null(data)) { + stop("Parameter 'data' cannot be NULL.") + } + if (is.null(quanti)) { + stop("Parameter 'quanti' is mandatory") + } + if (any(names(dim(data)) %in% 'sdate')) { + if (any(names(dim(data)) %in% 'ftime')) { + data <- MergeDims(data, merge_dims = c('ftime', 'sdate'), + rename_dim = 'time') + } + } + if (!(any(names(dim(data)) %in% 'time'))){ + stop("Parameter 'data' must have a temporal dimension named 'time'.") + } + if (any(names(dim(data)) %in% 'lat')) { + if (any(names(dim(data)) %in% 'lon')) { + data <- MergeDims(data, merge_dims = c('lon', 'lat'), + rename_dim = 'grid') + } + } + if (any(names(dim(data)) %in% 'latitude')) { + if (any(names(dim(data)) %in% 'longitude')) { + data <- MergeDims(data, merge_dims = c('longitude', 'latitude'), + rename_dim = 'grid') + } + } + if(!(any(names(dim(data)) %in% 'grid'))){ + stop("Parameter 'data' must have a spatial dimension named 'grid'.") + } + attractor <- Apply(data, target_dims = c('time', 'grid'), + fun = .proxiesattractor, + quanti = quanti , ncores = ncores) + # rename dimensions + attractor <- lapply(attractor, + FUN = function(x, dimname){ + names(dim(x))[dimname] <- 'time' + return(x)}, + dimname = which(names(dim(attractor[[1]])) == 'dim2')) + return(list(dim = attractor$dim, theta = attractor$theta)) +} + +.proxiesattractor <- function(data, quanti) { + # expected dimensions data: time and grid + logdista <- Apply(data, target_dims = 'grid', + fun = function(x, y){ + -log(colMeans((y - as.vector(x))^2))}, + y = t(data))[[1]] + + #Computation of theta + Theta <- function(logdista, quanti){ + #Compute the thheshold corresponding to the quantile + thresh <- quantile(logdista, quanti, na.rm = TRUE) + logdista[which(logdista == 'Inf')] <- NaN + Li <- which(as.vector(logdista) > as.numeric(thresh)) + #Length of each cluster + Ti <- diff(Li) + N <- length(Ti) + q <- 1 - quanti + Si <- Ti - 1 + Nc <- length(which(Si > 0)) + N <- length(Ti) + theta <- (sum(q * Si) + N + Nc - sqrt(((sum(q * Si) + N + Nc)^2) - + 8 * Nc * sum(q * Si))) / (2 * sum(q * Si)) + #Sort the exceedances + logdista <- sort(logdista) + #Find all the Peaks over Thresholds. + findidx <- which(as.vector(logdista) > as.numeric(thresh)) + if(length(findidx) < 1) { + stop("Parameter 'quanti' is too high for the length of the data provided.") + } + logextr <- logdista[findidx[[1]]:(length(logdista) - 1)] + #The inverse of the dimension is just the average of the exceedances + dim <- 1 /mean(as.numeric(logextr) - as.numeric(thresh)) + return(list(dim = dim, theta = theta)) + } + names(dim(logdista)) <- c('dim1', 'dim2') + proxies <- Apply(data = list(logdista = logdista), + target_dims = list('dim1'), fun = Theta, quanti = quanti) + + return(list(dim = proxies$dim, theta = proxies$theta)) +} + diff --git a/R/CST_SaveExp.R b/R/CST_SaveExp.R index 9c689ff7156d825d22cdb1bde8d4a6785c854dbc..f9290b18d739636b66fe61feeabb86e774e9b22e 100644 --- a/R/CST_SaveExp.R +++ b/R/CST_SaveExp.R @@ -13,6 +13,7 @@ #'folder tree: destination/experiment/variable/. By default the function #'creates and saves the data into the folder "CST_Data" in the working #'directory. +#'@param extra_string a character string to be include as part of the file name, for instance, to identify member or realization. It would be added to the file name between underscore characters. #' #'@seealso \code{\link{CST_Load}}, \code{\link{as.s2dv_cube}} and \code{\link{s2dv_cube}} #' @@ -29,7 +30,7 @@ #'} #' #'@export -CST_SaveExp <- function(data, destination = "./CST_Data") { +CST_SaveExp <- function(data, destination = "./CST_Data", extra_string = NULL) { if (!is.character(destination) & length(destination) > 1) { stop("Parameter 'destination' must be a character string of one element ", "indicating the name of the file (including the folder if needed) ", @@ -55,7 +56,8 @@ CST_SaveExp <- function(data, destination = "./CST_Data") { SaveExp(data = data$data, lon = data$lon, lat = data$lat, Dataset = names(data$Datasets), var_name = var_name, units = units, cdo_grid_name = cdo_grid_name, projection = projection, - startdates = sdates, Dates = time_values, destination) + startdates = sdates, Dates = time_values, destination, + extra_string = extra_string) } #'Save an experiment in a format compatible with CST_Load #'@description This function is created for compatibility with CST_Load/Load for saving post-processed datasets such as those calibrated of downscaled with CSTools functions @@ -73,8 +75,9 @@ CST_SaveExp <- function(data, destination = "./CST_Data") { #'@param cdo_grid_name a character string indicating the name of the grid e.g.: 'r360x181' #'@param projection a character string indicating the projection name #'@param destination a character string indicating the path where to store the NetCDF files +#'@param extra_string a character string to be include as part of the file name, for instance, to identify member or realization. #' -#'@return the function creates as many files as sdates per dataset. Each file could contain multiple members +#'@return the function creates as many files as sdates per dataset. Each file could contain multiple members. It would be added to the file name between underscore characters. #' The path will be created with the name of the variable and each Datasets. #' #'@import ncdf4 @@ -102,7 +105,8 @@ CST_SaveExp <- function(data, destination = "./CST_Data") { #'} #'@export SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates, - cdo_grid_name, projection, destination) { + cdo_grid_name, projection, destination, + extra_string = NULL) { dimname <- names(dim(data)) if (any(dimname == "ftime")) { dimname[which(dimname == "ftime")] <- "time" @@ -136,7 +140,11 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates, stop("Element 'data' in parameter 'data' has more than one 'sdate'", " dimension.") } - + if (!is.null(extra_string)) { + if (!is.character(extra_string)) { + stop("Parameter 'extra_string' must be a character string.") + } + } dataset_pos <- which(dimname == "dataset" | dimname == "dat") dims <- dim(data) if (length(dataset_pos) == 0) { @@ -227,7 +235,7 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates, NULL, 'time'), fun = .saveExp, var_name = var_name, units = units, dims_var = dims_var, cdo_grid_name = cdo_grid_name, projection = projection, - destination = path) + destination = path, extra_string = extra_string) } } @@ -252,7 +260,7 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates, #Dates <- as.Date(c("1900-11-01", "1900-12-01", "1901-01-01", "1901-02-01", "1901-03-01")) #.saveExp(data, sdate, Dates, var_name, units, dims_var, cdo_grid_name = 'r360x181', projection = 'none', destination) .saveExp <- function(data, sdate, Dates, var_name, units, dims_var, - cdo_grid_name, projection, destination) { + cdo_grid_name, projection, destination, extra_string) { dim_names <- names(dim(data)) if (any(dim_names != c('longitude', 'latitude', 'member', 'time'))) { data <- Reorder(data, c('longitude', 'latitude', 'member', 'time')) @@ -266,7 +274,11 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates, datanc <- ncvar_def(name = var_name, units = units, dim = dims_var, missval = -99999) - file_name <- paste0(var_name, "_", sdate, ".nc") + if (is.null(extra_string)) { + file_name <- paste0(var_name, "_", sdate, ".nc") + } else { + file_name <- paste0(var_name, "_", extra_string, "_", sdate, ".nc") + } full_filename <- file.path(destination, file_name) file_nc <- nc_create(full_filename, datanc) ncvar_put(file_nc, datanc, data) diff --git a/R/PlotCombinedMap.R b/R/PlotCombinedMap.R index 2aee71b8ef7a3861a3834e05f0110ea10bc1751d..8af4a14db9f69e258e1eb5bf28961cace1be9def 100644 --- a/R/PlotCombinedMap.R +++ b/R/PlotCombinedMap.R @@ -14,8 +14,15 @@ #'@param col_unknown_map Colour to use to paint the grid cells for which a map is not possible to be chosen according to 'map_select_fun' or for those values that go beyond 'display_range'. Takes the value 'white' by default. #'@param mask Optional numeric array with dimensions (latitude, longitude), with values in the range [0, 1], indicating the opacity of the mask over each grid point. Cells with a 0 will result in no mask, whereas cells with a 1 will result in a totally opaque superimposed pixel coloured in 'col_mask'. #'@param col_mask Colour to be used for the superimposed mask (if specified in 'mask'). Takes the value 'grey' by default. +#'@param dots Array of same dimensions as 'var' or with dimensions +#' c(n, dim(var)), where n is the number of dot/symbol layers to add to the +#' plot. A value of TRUE at a grid cell will draw a dot/symbol on the +#' corresponding square of the plot. By default all layers provided in 'dots' +#' are plotted with dots, but a symbol can be specified for each of the +#' layers via the parameter 'dot_symbol'. #'@param bar_titles Optional vector of character strings providing the titles to be shown on top of each of the colour bars. #'@param legend_scale Scale factor for the size of the colour bar labels. Takes 1 by default. +#'@param cex_bar_titles Scale factor for the sizes of the bar titles. Takes 1.5 by default. #'@param fileout File where to save the plot. If not specified (default) a graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, bmp and tiff #'@param width File width, in the units specified in the parameter size_units (inches by default). Takes 8 by default. #'@param height File height, in the units specified in the parameter size_units (inches by default). Takes 5 by default. @@ -61,7 +68,9 @@ PlotCombinedMap <- function(maps, lon, lat, brks = NULL, cols = NULL, col_unknown_map = 'white', mask = NULL, col_mask = 'grey', + dots = NULL, bar_titles = NULL, legend_scale = 1, + cex_bar_titles = 1.5, fileout = NULL, width = 8, height = 5, size_units = 'in', res = 100, ...) { @@ -280,7 +289,17 @@ PlotCombinedMap <- function(maps, lon, lat, stop("Parameter 'mask' must have dimensions c(lat, lon).") } } - + # Check dots + if (!is.null(dots)) { + if (length(dim(dots)) != 2) { + stop("Parameter 'mask' must have two dimensions.") + } + if ((dim(dots)[1] != dim(maps)[lat_dim]) || + (dim(dots)[2] != dim(maps)[lon_dim])) { + stop("Parameter 'mask' must have dimensions c(lat, lon).") + } + } + #---------------------- # Identify the most likely map #---------------------- @@ -327,6 +346,9 @@ PlotCombinedMap <- function(maps, lon, lat, if (!is.null(mask)){ mask <- mask[nlat:1, ] } + if (!is.null(dots)){ + dots <- dots[nlat:1,] + } } #---------------------- @@ -353,7 +375,7 @@ PlotCombinedMap <- function(maps, lon, lat, tbrks <- c(-1, brks_norm + rep(1:nmap, each = length(brks))) PlotEquiMap(var = ml_map, lon = lon, lat = lat, brks = tbrks, cols = tcols, drawleg = FALSE, - filled.continents = FALSE, ...) + filled.continents = FALSE, dots = dots, ...) #---------------------- # Add overplot on top @@ -402,7 +424,7 @@ PlotCombinedMap <- function(maps, lon, lat, draw_separators = TRUE, extra_margin = c(2, 0, 2, 0), label_scale = legend_scale * 1.5) if (!is.null(bar_titles)) { - mtext(bar_titles[[k]], 3, line = -3, cex = 1.5) + mtext(bar_titles[[k]], 3, line = -3, cex = cex_bar_titles) } } diff --git a/R/Predictability.R b/R/Predictability.R new file mode 100644 index 0000000000000000000000000000000000000000..12cd6b410f26b197da67df65e82f42f626d29668 --- /dev/null +++ b/R/Predictability.R @@ -0,0 +1,163 @@ +#'@rdname Predictability +#'@title Computing scores of predictability using two dynamical proxies +#'based on dynamical systems theory. +#' +#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} +#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} +#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it} +#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +#' +#'@description This function divides in terciles the two dynamical proxies +#'computed with CST_ProxiesAttractor or ProxiesAttractor. These terciles will +#'be used to measure the predictability of the system in dyn_scores. When the +#'local dimension 'dim' is small and the inverse of persistence 'theta' is +#'small the predictability is high, and viceversa. +#' +#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., +#'and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large +#'scale atmospheric predictability.Nature Communications, 10(1), 1316. +#'DOI = https://doi.org/10.1038/s41467-019-09305-8 " +#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +#' Dynamical proxies of North Atlantic predictability and extremes. +#' Scientific Reports, 7-41278, 2017. +#' +#'@param dim An array of N named dimensions containing the local dimension as +#'the output of CST_ProxiesAttractor or ProxiesAttractor. +# +#'@param theta An array of N named dimensions containing the inverse of the +#'persistence 'theta' as the output of CST_ProxiesAttractor or ProxiesAttractor. +#' +#'@param ncores The number of cores to use in parallel computation +#' +#'@return A list of length 2: +#' \itemize{ +#' \item\code{pred.dim} {a list of two lists 'qdim' and 'pos.d'. The 'qdim' list +#'contains values of local dimension 'dim' divided by terciles: +#'d1: lower tercile (high predictability), +#'d2: middle tercile, +#'d3: higher tercile (low predictability) +#'The 'pos.d' list contains the position of each tercile in parameter 'dim'} +#' +#' \item\code{pred.theta} {a list of two lists 'qtheta' and 'pos.t'. +#'The 'qtheta' list contains values of the inverse of persistence 'theta' +#'divided by terciles: +#'th1: lower tercile (high predictability), +#'th2: middle tercile, +#'th3: higher tercile (low predictability) +#'The 'pos.t' list contains the position of each tercile in parameter 'theta'} +#'} +#'@return dyn_scores values from 0 to 1. A dyn_score of 1 indicates the highest +#'predictability. +#' +#'@examples +#'# Creating an example of matrix dat(time,grids): +#'m <- matrix(rnorm(2000) * 10, nrow = 50, ncol = 40) +#'names(dim(m)) <- c('time', 'grid') +#'# imposing a threshold +#' quanti <- 0.90 +#'# computing dyn_scores from parameters dim and theta of the attractor +#' attractor <- ProxiesAttractor(dat = m, quanti = 0.60) +#' predyn <- Predictability(dim = attractor$dim, theta = attractor$theta) +#'@export +#' +Predictability<- function(dim, theta, ncores = NULL) { + # if (!inherits(dim, 's2dv_cube')) { + # stop("Parameter 'dim' must be of the class 's2dv_cube', ", + # "as output by CSTools::CST_Load.") + # } + # if (!inherits(theta, 's2dv_cube')) { + # stop("Parameter 'theta' must be of the class 's2dv_cube', ", + # "as output by CSTools::CST_Load.") + # } + if (any(names(dim(dim)) %in% 'sdate')) { + if (any(names(dim(dim)) %in% 'ftime')) { + dim <- MergeDims(dim, merge_dims = c('ftime', 'sdate'), + rename_dim = 'time') + } + } + if (!(any(names(dim(dim)) %in% 'time'))){ + stop("Parameter 'dim' must have a temporal dimension named 'time'.") + } + + if (any(names(dim(theta)) %in% 'sdate')) { + if (any(names(dim(theta)) %in% 'ftime')) { + theta <- MergeDims(theta, merge_dims = c('ftime', 'sdate'), + rename_dim = 'time') + } + } + if (!(any(names(dim(theta)) %in% 'time'))){ + stop("Parameter 'data' must have a temporal dimension named 'time'.") + } + + pred <- Apply(list(dim, theta), target_dims = 'time', + fun = .predictability, + ncores = ncores) + dim(pred$dyn_scores) <- dim(theta) + return(list(pred.dim = list(qdim = list(pred$qdim.d1,pred$qdim.d2,pred$qdim.d3), + pos.d = list(pred$pos.d1,pred$pos.d2,pred$pos.d3)), + pred.theta = list(qtheta = list(pred$qtheta.th1,pred$qtheta.th2,pred$qtheta.th3), + pos.t = list(pred$pos.th1,pred$pos.th2,pred$pos.th3)), + dyn_scores = pred$dyn_scores)) +} +.predictability <- function(dim, theta) { + if (is.null(dim)) { + stop("Parameter 'dim' cannot be NULL.") + } + if (is.null(theta)) { + stop("Parameter 'theta' cannot be NULL.") + } + if (length(dim) != length(theta)) { + stop("Parameters 'dim' and 'theta' must have the same length") + } + pos <- c(1:length(dim)) + + # dim + qd1 <- quantile(dim, 0.33, na.rm = TRUE) + qd3 <- quantile(dim, 0.66, na.rm = TRUE) + d3 <- which(dim >= qd3) + d1 <- which(dim <= qd1) + d2 <- which(dim > qd1 & dim < qd3) + qdim <- list(d1 = dim[d1], d2 = dim[d2], d3 = dim[d3]) + pos.d <- list(d1=d1,d2=d2,d3=d3) + + #theta + qt1 <- quantile(theta, 0.33, na.rm = TRUE) + qt3 <- quantile(theta, 0.66, na.rm = TRUE) + th3 <- which(theta >= qt3) + th1 <- which(theta <= qt1) + th2 <- which(theta > qt1 & theta < qt3) + qtheta <- list(th1 = theta[th1], th2 = theta[th2], th3 = theta[th3]) + pos.t <- list(th1 = th1, th2 = th2, th3 = th3) + + #scores position + d1th1 <- pos.d$d1[which(pos.d$d1 %in% pos.t$th1)] + d1th2 <- pos.d$d1[which(pos.d$d1 %in% pos.t$th2)] + d2th1 <- pos.d$d2[which(pos.d$d2 %in% pos.t$th1)] + d1th3 <- pos.d$d1[which(pos.d$d1 %in% pos.t$th3)] + d3th1 <- pos.d$d3[which(pos.d$d3 %in% pos.t$th1)] + d2th2 <- pos.d$d2[which(pos.d$d2 %in% pos.t$th2)] + d2th3 <- pos.d$d2[which(pos.d$d2 %in% pos.t$th3)] + d3th2 <- pos.d$d3[which(pos.d$d3 %in% pos.t$th2)] + d3th3 <- pos.d$d3[which(pos.d$d3 %in% pos.t$th3)] + #scores values + dyn_scores <- c(1:length(dim)) + dyn_scores[which(pos %in% d1th1)]<- 9/9 + dyn_scores[which(pos %in% d1th2)]<- 8/9 + dyn_scores[which(pos %in% d2th1)]<- 7/9 + dyn_scores[which(pos %in% d1th3)]<- 6/9 + dyn_scores[which(pos %in% d3th1)]<- 5/9 + dyn_scores[which(pos %in% d2th2)]<- 4/9 + dyn_scores[which(pos %in% d2th3)]<- 3/9 + dyn_scores[which(pos %in% d3th2)]<- 2/9 + dyn_scores[which(pos %in% d3th3)]<- 1/9 + +return(list(qdim.d1 = dim[d1], qdim.d2 = dim[d2], qdim.d3 = dim[d3], + pos.d1 = d1, pos.d2 = d2, pos.d3 =d3, + qtheta.th1 = theta[th1], qtheta.th2 = theta[th2], qtheta.th3 = theta[th3], + pos.th1 = th1, pos.th2 = th2, pos.th3 = th3, + dyn_scores = dyn_scores)) +} + + + + diff --git a/inst/doc/UseCase1_WindEvent_March2018.R b/inst/doc/UseCase1_WindEvent_March2018.R new file mode 100644 index 0000000000000000000000000000000000000000..d3bbb9363ba0c9058efa573bae6f0a9512163a80 --- /dev/null +++ b/inst/doc/UseCase1_WindEvent_March2018.R @@ -0,0 +1,211 @@ +#!/usr/bin/env Rscript +rm(list=ls()); gc(); + +### Creation date: 3rd June 2021 +# Author: N. Pérez-Zanón +# Refer CSTools package and manuscript when using it. +# ---------------------------------------- +# Wind speed bias adjustment for assessment of an extreme event +# The System5-ECMWF downloaded from C3S seasonal forecasts initialized +# in December 2017, January 2018 and February 2018 +# This code includes the bias adjustent and the results visualization +# ---------------------------------------- + +#library(CSTools) +library(s2dv) +library(ragg) +library(multiApply) +output_dir <- "/esarchive/scratch/nperez/CSTools_manuscript/v20210603/" + + +exp_path <- list(name = "ECMWFS5", + path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_f6h/$VAR_NAME$_$START_DATE$.nc") +obs_path <- list(name = "ERA5", + path = "/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h/$VAR_NAME$_$YEAR$$MONTH$.nc") +#source("/esarchive/scratch/nperez/git/cstools/R/CST_BiasCorrection.R") +#library(multiApply) + # Target months March (3) + # Assess forecast from 1 to 3 months in advance +months_in_advance <- c('02', '01', '12') +wind_fsct_BC <- list() +wind_ref_terciles <- NULL +wind_ref_extremes <- NULL +wind_thres_latlon <- NULL + # Observations March 2018 + wind_obs <- CSTools::CST_Load(var = 'windagl100', obs = list(obs_path), + sdates = '20180301', nmember = 1, + leadtimemin = 1, leadtimemax = 1, + storefreq = "monthly", sampleperiod = 1, + latmin = 36, latmax = 44, lonmin = -10, lonmax = 4, +# latmin = 42, latmax = 44, lonmin = -10, lonmax = 1, + output = 'lonlat', nprocs = 1, grid = 'r360x181') + # For each month in advance: +for (mm in 1:3) { + print("Three initializations:") + print(mm) + print(paste('Initialization', months_in_advance[mm])) + # Generate the start dates of hindcast period + year <- ifelse(mm == 3, 2017, 2018) + print(paste('Hindcast period until', year - 1)) + hcst_sdates <- paste0(1993:(year - 1), months_in_advance[mm], '01') + wind_hcst <- CSTools::CST_Load(var = 'sfcWind', exp = list(exp_path), + sdates = hcst_sdates, nmember = 25, + leadtimemin = mm + 1, leadtimemax = mm + 1, + storefreq = "monthly", sampleperiod = 1, + latmin = 36, latmax = 44, lonmin = -10, lonmax = 4, + # latmin = 42, latmax = 44, lonmin = -10, lonmax = 1, + output = 'lonlat', nprocs = 1) + str(wind_hcst$Dates) + dim(wind_hcst$data) + fcst_sdates <- paste0(year, months_in_advance[mm], '01') + wind_fcst <- CSTools::CST_Load(var = 'sfcWind', exp = list(exp_path), + sdates = fcst_sdates, nmember = 25, + leadtimemin = mm + 1, leadtimemax = mm + 1, + storefreq = "monthly", sampleperiod = 1, + latmin = 36, latmax = 44, lonmin = -10, lonmax = 4, + # latmin = 42, latmax = 44, lonmin = -10, lonmax = 1, + + output = 'lonlat', nprocs = 1) + str(wind_fcst$Dates) + dim(wind_fcst$data) + + wind_ref <- CSTools::CST_Load(var = 'windagl100', obs = list(obs_path), + sdates = hcst_sdates, nmember = 1, + leadtimemin = mm + 1, leadtimemax = mm + 1, + storefreq = "monthly", sampleperiod = 1, + latmin = 36, latmax = 44, lonmin = -10, lonmax = 4, +# latmin = 42, latmax = 44, lonmin = -10, lonmax = 1, + output = 'lonlat', nprocs = 1, + grid = 'r360x181') + str(wind_ref$Dates) + dim(wind_ref$data) + print(wind_ref$Dates$start) + + wind_ref_terciles <- rbind(wind_ref_terciles, + quantile(MeanDims(wind_ref$data, c('lat', 'lon')), c(0.3, 0.6))) + wind_ref_extremes <- rbind(wind_ref_extremes, + quantile(MeanDims(wind_ref$data, c('lat', 'lon')), c(0.1, 0.9))) + source("/esarchive/scratch/nperez/git/cstools/R/CST_BiasCorrection.R") + wind_fsct <- CST_BiasCorrection(exp = wind_hcst, + obs = wind_ref, + exp_cor = wind_fcst) + wind_fsct_BC[[mm]] <- wind_fsct + # ----------------------------------------------- + # PLOTTING: + # PlotMostlikely + thres <- drop(Apply(list(wind_ref$data), target_dims = 'sdate', fun = function(x) { + quantile(x, c(1/3, 2/3))}, output_dims = 'probs')$output) + PB <- Apply(list(wind_fsct$data, thres), target_dims = list('member', 'probs'), + fun = function(x, z) { + result <- unlist(lapply(x, function(y) { + if (y <= z[1]) { + res <- c(1, 0, 0) + } else if (y <= z[2]) { + res <- c(0, 1, 0) + } else { + res <- c(0, 0, 1) + } + return(res)})) + dim(result) <- c(bin = 3, member = 25) + return(result) + })$output1 + Mean_PB <- drop(MeanDims(PB, 'member')) + + observed_tercile <- Apply(list(wind_obs$data, thres, Mean_PB), + target_dims = list(NULL, 'probs', 'bin'), + fun = function(x, y, z) { + val <- which.max(z) + if (val == 3) { + dot <- ifelse(x >= y[2], 1, 0) + } else if (val == 2) { + dot <- ifelse(x < y[2] && x >= y[1], 1, 0) + } else if (val == 1) { + dot <- ifelse(x < y[1], 1, 0) + } else { + stop('what') + } + return(dot) + })$output1 + + filtered_obs_terciles <- Apply(list(Mean_PB, observed_tercile), + target_dims = list('bin', NULL), + fun = function(x,y) { + if (sum(duplicated(x)) == 1) { + y <- 0 + } + return(y) })$output1 + + wind_thres_latlon <- abind::abind(wind_thres_latlon, thres, along = 4) + source("/esarchive/scratch/nperez/git/cstools/R/PlotCombinedMap.R") + source("/esarchive/scratch/nperez/git/cstools/R/PlotMostLikelyQuantileMap.R") + agg_png(paste0(output_dir, "Wind_MostLikely_", mm, "_obstercile.png"), + width = 1050, height = 1000, units = 'px', res = 144) + PlotMostLikelyQuantileMap(probs = Mean_PB, lon = wind_fsct$lon, + lat = wind_fsct$lat, sizetit = 1.5, + intylat = 2, intxlon = 2, + coast_width = 1.5, legend_scale = 0.8, + cat_dim = 'bin', dot_size = 2.5, + axes_label_scale = 1.6, degree_sym = TRUE, + dots = filtered_obs_terciles[,,1,1,1,1], + toptitle = c(paste0("Initialized on ", + month.name[as.numeric(months_in_advance[mm])]))) + dev.off() +} + +visual <- data.frame(dec = as.vector(MeanDims(wind_fsct_BC[[3]]$data, c('lat', 'lon'))), + jan = as.vector(MeanDims(wind_fsct_BC[[2]]$data, c('lat', 'lon'))), + feb = as.vector(MeanDims(wind_fsct_BC[[1]]$data, c('lat', 'lon')))) + + wind_obs_areave <- CSTools::CST_Load(var = 'windagl100', obs = list(obs_path), + sdates = '20180301', nmember = 1, + leadtimemin = 1, leadtimemax = 1, + storefreq = "monthly", sampleperiod = 1, + latmin = 36, latmax = 44, lonmin = -10, lonmax = 4, +# latmin = 42, latmax = 44, lonmin = -10, lonmax = 1, + output = 'areave', nprocs = 1) + +print("IS DATA LOADED") + +print("Wait") +agg_png(paste0(output_dir, "Wind_PlotForecast_IP.png"), + width = 1000, height = 1000, units = 'px',res = 150) +CSTools::PlotForecastPDF(visual, tercile.limits = wind_ref_terciles, + extreme.limits = wind_ref_extremes, + var.name = "Wind Speed 100 m (m/s)", + title = "Bias Corrected forecasts at IP", + fcst.names = c("December", "January", "February"), + obs = as.vector(wind_obs_areave$data)) +dev.off() + +# Plotting observed terciles: +names(dim(wind_thres_latlon)) <- c('thres', 'lat', 'lon', 'sdate') +wind_thres_latlon <- ClimProjDiags::Subset(wind_thres_latlon, indices = 1, along = 'sdate') +wind_obs_obstercile <- Apply(list(wind_obs$data, wind_thres_latlon), + target_dims = list(NULL, 'thres'), + fun = function(x, tercile) { + if (x <= tercile[1]) { + res <- 1 + } else if (x > tercile[2]) { + res <- 3 + } else { + res <- 2 + } + return(res) + })$output1 + +agg_png(paste0(output_dir, "MostLikely_Observed_obstercile.png"), + width = 1000, height = 1000, units = 'px', res = 144) + +s2dv::PlotEquiMap(wind_obs_obstercile, + lon = wind_obs$lon, lat = wind_obs$lat, + brks = c(0,1,2,3), + cols = c("#6BAED6FF", "#FFEDA0FF", "#FC4E2AFF"), + intylat = 2, intxlon = 2, + coast_width = 1.5, filled.continents = FALSE, + toptitle = c("Observed terciles March 2018")) +dev.off() + +# All gridpoints are above normal observed tercile. + +print("DONE") + diff --git a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF100.R b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF100.R new file mode 100644 index 0000000000000000000000000000000000000000..146382dc8983caa3873a8b13665214b6d0dfb7d0 --- /dev/null +++ b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF100.R @@ -0,0 +1,160 @@ +#!/usr/bin/env Rscript +rm(list=ls()); gc(); + +### Creation date: 9th June 2021 +# Author: N. Pérez-Zanón +# Adapted from the original version Terzago et al. 2020 +# https://drive.google.com/file/d/1qp2gbtKdBl4XmsyOeaEhFENwpeUuJwkf/view +# Refer CSTools package and manuscript when using it. +# ---------------------------------------- +# This code is using divided in 4 steps plus visualization of the results +# Taking advantage of +# STEP 1: Compute Slope values for RainFARM downscaling method +# STEP 2: Apply Quantile Mapping Correction to Seasonal Precipitation Forecast +# STEP 3: Compute Weights for RainFARM downscaling method +# STEP 4: Apply RainFARM downscaling method +# ---------------------------------------- +# Note: STEP 3 requires a machine with internet connexion. +# In this file, the lines are commented since they have been run and the +# result saved on disk, then the result is loaded. +# ---------------------------------------- +# +# Load required libraries and setup output directory: +library(CSTools) +library(ClimProjDiags) +library(zeallot) +library(ragg) +dir_output <- '/esarchive/scratch/nperez/CSTools_manuscript/v20210603/' #slash end + +# -------------------------------------------- +# STEP 1: +# -------------------------------------------- +era5 <- list(name = 'era5', + path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc') +years <- unlist(lapply(1993:2018, function(x){ + paste0(x, sprintf("%02d",1:12), '01')})) +era5 <- CST_Load(var = 'prlr', + exp = list(era5), + sdates = years, nmember = 1, + storefreq = "daily", sampleperiod = 1, + latmin = 37.5, latmax = 53.25, lonmin = 2.5, lonmax = 18.25, + output = 'lonlat', nprocs = 1) +era5$data <- era5$data * 24 * 3600 * 1000 # with or without this line -> same result +era5 <- CST_SplitDim(era5, split_dim = 'sdate', indices = rep(1:12, 26)) +slope <- CST_RFSlope(era5, time_dim = c('sdate', 'ftime'), kmin = 5) +save(slope, file = paste0(dir_output, 'Slope.RDS'), version = 2) +slope_plot <- slope + +# -------------------------------------------- +# STEP 2: +# -------------------------------------------- +StartDates <- paste0(1993:2018, '1101') +exp <- list(name = 'ecmwfS5', + path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc") +obs <- list(name = 'era5', + path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc') +#obs <- list(name = 'era5', path = '/esarchive/scratch/nperez/ERA5/daily_total_prlr_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc') + +c(exp, obs) %<-% CST_Load(var = 'prlr', exp = list(exp), obs = list(obs), + sdates = StartDates, nmember = 25, + storefreq = "daily", sampleperiod = 1, + latmin = 42, latmax = 49, lonmin = 4, lonmax = 11, + output = 'lonlat', nprocs = 1) +exp <- CST_SplitDim(exp, split_dim = c('ftime')) +obs <- CST_SplitDim(obs, split_dim = c('ftime')) +exp$data <- exp$data * 24 * 3600 * 1000 +obs$data <- obs$data * 24 * 3600 * 1000 +exp$data[which(exp$data < 0)] <- 0 +exp.qm <- CST_QuantileMapping(exp, obs, method = "QUANT", + wet.day = FALSE, + sample_dims = c('member', 'sdate', 'ftime'), + ncores = 4) +save(exp.qm, file = paste0(dir_output, 'ExpQM.RDS'), version = 2) + +# -------------------------------------------- +# STEP 3: +# -------------------------------------------- +#library(raster); library(s2dv); library(CSTools) +#worldclim <- getData("worldclim", var = "prec", res = 0.5, lon = 5, lat = 45) +#wc_month <- lapply(1:12, FUN = function(x) { +# res <- crop(worldclim[[x]], +# extent(3.5, 11.5, 41.5, 49.5)) +# res <- as.array(res) +# names(dim(res)) <- c('lat', 'lon', 'month') +# return(res) +# }) +#xy <- xyFromCell(crop(worldclim[[1]], extent(3.5, 11.5, 41.5, 49.5)), +# 1:length(crop(worldclim[[1]], extent(3.5, 11.5, 41.5, 49.5)))) +#lons <- unique(xy[,1]) +#lats <- unique(xy[,2]) +#wc_month <- unlist(wc_month) +#dim(wc_month) <- c(lat = length(lats), lon = length(lons), monthly = 12) +#wc_month <- Reorder(wc_month, c('lon', 'lat', 'monthly')) +#wc_month <- s2dv_cube(data = wc_month, lon = lons, lat = lats, +# Datasets = 'WorldClim') +#weight <- CST_RFWeights(wc_month, lon = 4:11, lat = 49:42, nf = 100) +#save(weight, file = paste0(dir_output, 'weightsRF100.RDS')) + +load(paste0(dir_output, 'weightsRF100.RDS')) + +# -------------------------------------------- +# Visualization +# -------------------------------------------- +agg_png(paste0(dir_output, "RF100_WeightsDec.png"), + width = 1000, height = 1100, units = 'px',res = 144) +PlotEquiMap(weight$data[,,12], lon = weight$lon, lat = weight$lat, + filled.continents = FALSE, title_scale = 1, + intylat = 2, intxlon = 2, + toptitle = 'December Weights RF 100') +dev.off() + +# -------------------------------------------- +# STEP 4: +# -------------------------------------------- + +weights <- Subset(weight$data, along = 'monthly', indices = c(11, 12, 1:6)) +slope <- Subset(slope, along = 'monthly', indices = c(11, 12, 1:6), + drop = 'non-selected') +k = 1 # To create the member ID +for (realizations in 1:10) { + for (member in 1:25) { + result <- exp.qm # to store the data + result$data <- NULL + for (month in 1:8) { + data <- exp.qm # to take the correct piece of data + data$data <- data$data[1, member, , , , , month] + fs <- CST_RainFARM(data, nf = 100, + weights = weights, slope = slope[month], + kmin = 1, nens = 1, verbose = TRUE, + nprocs = 8, + drop_realization = TRUE) + result$data <- abind::abind(result$data, fs$data, along = 5) + if (month == 2 & member == 1 & realization == 1) { + # ---------------------------- + # Visualization: + # ---------------------------- + agg_png(paste0(dir_output, "RF100_Down_11dec.png"), + width = 1000, height = 1100, units = 'px',res = 144) + PlotEquiMap(fs$data[1,11,,],lon = fs$lon, lat = fs$lat, + filled.continents = FALSE, bar_limits = c(0,40), + intylat = 2, intxlon = 2, title_scale = 1, + triangle_ends = c(TRUE, FALSE), + toptitle = 'Downsacaled RF 100', units = 'precipitation (mm)') + dev.off() + } + result$lon <- fs$lon + result$lat <- fs$lat + result <- CST_MergeDims(result, merge_dims = c("ftime", "monthly"), + na.rm = TRUE) + result$Dataset <- paste0('RF100_ECMWFC3S_QM_member_', member, '_real_', + realizations) + result$Dates[[1]] <- exp$Dates[[1]] + CST_SaveExp(result, destination = dir_output, + extra_string = paste0('member', k)) + gc() + k = k + 1 + rm(list= list('fs', 'result', 'data')) + } +} + + diff --git a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R new file mode 100644 index 0000000000000000000000000000000000000000..ee2df1943eef5547d21401f55c8fe10d2853132d --- /dev/null +++ b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R @@ -0,0 +1,215 @@ +#!/usr/bin/env Rscript +rm(list=ls()); gc(); + +### Creation date: 3rd June 2021 +# Author: N. Pérez-Zanón +# Adapted from the original version Terzago et al. 2020 +# https://drive.google.com/file/d/1qp2gbtKdBl4XmsyOeaEhFENwpeUuJwkf/view +# Refer CSTools package and manuscript when using it. +# ---------------------------------------- +# This code is divided in 4 steps plus visualization of the results +# STEP 1: Compute Slope values for RainFARM downscaling method +# STEP 2: Apply Quantile Mapping Correction to Seasonal Precipitation Forecast +# STEP 3: Compute Weights for RainFARM downscaling method +# STEP 4: Apply RainFARM downscaling method +# ---------------------------------------- +# Note: STEP 3 requires a machine with internet connexion. +# In this file, the lines are commented since they have been run and the +# result saved on disk, then the result is loaded. +# ---------------------------------------- +# +# Load required libraries and setup output directory: +library(CSTools) +library(ClimProjDiags) +library(zeallot) +library(ragg) +dir_output <- '/esarchive/scratch/nperez/CSTools_manuscript/v20210603/' #slash end +# -------------------------------------------- +# STEP 1: +# -------------------------------------------- +era5 <- list(name = 'era5', + path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc') +years <- unlist(lapply(1993:2018, function(x){ + paste0(x, sprintf("%02d",1:12), '01')})) +era5 <- CST_Load(var = 'prlr', + exp = list(era5), + sdates = years, nmember = 1, + storefreq = "daily", sampleperiod = 1, + latmin = 37.5, latmax = 53.25, lonmin = 2.5, lonmax = 18.25, + output = 'lonlat', nprocs = 1) +era5$data <- era5$data * 24 * 3600 * 1000 # with or without this line -> same result +era5 <- CST_SplitDim(era5, split_dim = 'sdate', indices = rep(1:12, 26)) +slope <- CST_RFSlope(era5, time_dim = c('sdate', 'ftime'), kmin = 5) +save(slope, file = paste0(dir_output, 'Slope.RDS'), version = 2) +slope_plot <- slope + +# -------------------------------------------- +# STEP 2: +# -------------------------------------------- +StartDates <- paste0(1993:2018, '1101') +exp <- list(name = 'ecmwfS5', + path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc") +obs <- list(name = 'era5', + path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc') +#obs <- list(name = 'era5', path = '/esarchive/scratch/nperez/ERA5/daily_total_prlr_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc') + +c(exp, obs) %<-% CST_Load(var = 'prlr', exp = list(exp), obs = list(obs), + sdates = StartDates, nmember = 25, + storefreq = "daily", sampleperiod = 1, + latmin = 42, latmax = 49, lonmin = 4, lonmax = 11, + output = 'lonlat', nprocs = 1) +exp <- CST_SplitDim(exp, split_dim = c('ftime')) +obs <- CST_SplitDim(obs, split_dim = c('ftime')) +exp$data <- exp$data * 24 * 3600 * 1000 +obs$data <- obs$data * 24 * 3600 * 1000 +exp$data[which(exp$data < 0)] <- 0 +exp.qm <- CST_QuantileMapping(exp, obs, method = "QUANT", + wet.day = FALSE, + sample_dims = c('member', 'sdate', 'ftime'), + ncores = 4) +save(exp.qm, file = paste0(dir_output, 'ExpQM.RDS'), version = 2) +save(exp, file = paste0(dir_output, 'Exp.RDS'), version = 2) +# -------------------------------------------- +# STEP 3: +# -------------------------------------------- +#library(raster); library(s2dv); library(CSTools) +#worldclim <- getData("worldclim", var = "prec", res = 0.5, lon = 5, lat = 45) +#wc_month <- lapply(1:12, FUN = function(x) { +# res <- crop(worldclim[[x]], +# extent(3.5, 11.5, 41.5, 49.5)) +# res <- as.array(res) +# names(dim(res)) <- c('lat', 'lon', 'month') +# return(res) +# }) +#xy <- xyFromCell(crop(worldclim[[1]], extent(3.5, 11.5, 41.5, 49.5)), +# 1:length(crop(worldclim[[1]], extent(3.5, 11.5, 41.5, 49.5)))) +#lons <- unique(xy[,1]) +#lats <- unique(xy[,2]) +#wc_month <- unlist(wc_month) +#dim(wc_month) <- c(lat = length(lats), lon = length(lons), monthly = 12) +#wc_month <- Reorder(wc_month, c('lon', 'lat', 'monthly')) +#wc_month <- s2dv_cube(data = wc_month, lon = lons, lat = lats, +# Datasets = 'WorldClim') +#weight <- CST_RFWeights(wc_month, lon = 4:11, lat = 49:42, nf = 4) +#save(weight, file = paste0(dir_output, 'weightsRF4.RDS')) + +load(paste0(dir_output, 'weightsRF4.RDS')) + +# -------------------------------------------- +# STEP 4: +# -------------------------------------------- +Rprof() +weights <- Subset(weight$data, along = 'monthly', indices = c(11, 12, 1:6)) +slope <- Subset(slope, along = 'monthly', indices = c(11, 12, 1:6), + drop = 'non-selected') +fs <- CST_RainFARM(exp.qm, nf = 4, + weights = weights, slope = slope, + kmin = 1, nens = 10, verbose = TRUE, + time_dim = c("member", "sdate", "ftime"), nprocs = 8, + drop_realization = TRUE) +newfs <- CST_MergeDims(fs, merge_dims = c("ftime", "monthly"), + na.rm = TRUE) + +newfs$Dates[[1]] <- exp$Dates[[1]] +CST_SaveExp(newfs, destination = paste0(dir_output, 'RF4/')) + +Rprof(NULL) +profile.info <- summaryRprof(paste0(dir_output, "Rprof.out")) +# -------------------------------------------- +# Visualization +# -------------------------------------------- +library(s2dv) +agg_png(paste0(dir_output, "EXP_11dec.png"), + width = 800, height = 900, units = 'px',res = 144) +PlotEquiMap(exp$data[1,1,1,11,,,2],lon = exp$lon, lat = exp$lat, + filled.continents = FALSE, bar_limits = c(0,40), + intylat = 2, intxlon = 2, title_scale = 0.8, + bar_label_scale = 1.3, axes_label_scale = 1.2, + triangle_ends = c(TRUE, FALSE), degree_sym = TRUE, units_scale = 1.5, + toptitle = 'SEAS5', units = 'precipitation (mm)') +dev.off() +agg_png(paste0(dir_output, "EXPQM_11dec.png"), + width = 800, height = 900, units = 'px',res = 144) +PlotEquiMap(exp.qm$data[1,1,1,11,,,2],lon = exp$lon, lat = exp$lat, + filled.continents = FALSE, bar_limits = c(0,40), + intylat = 2, intxlon = 2, title_scale = 0.8, + bar_label_scale = 1.3, axes_label_scale = 1.2, + triangle_ends = c(TRUE, FALSE), degree_sym = TRUE, units_scale = 1.5, + toptitle = 'Bias Adjusted', units = 'precipitation (mm)') +dev.off() +agg_png(paste0(dir_output, "RF4_Down_11dec.png"), + width = 800, height = 900, units = 'px',res = 144) +PlotEquiMap(fs$data[1,1,1,11,,,2],lon = fs$lon, lat = fs$lat, + filled.continents = FALSE, bar_limits = c(0,40), + intylat = 2, intxlon = 2, title_scale = 0.8, + bar_label_scale = 1.3, axes_label_scale = 1.2, + triangle_ends = c(TRUE, TRUE), degree_sym = TRUE, units_scale = 1.5, + toptitle = 'Downscaled nf 4', units = 'precipitation (mm)') +dev.off() +agg_png(paste0(dir_output, "RF4_WeightsDec.png"), + width = 800, height = 900, units = 'px',res = 144) +PlotEquiMap(weight$data[,,12], lon = weight$lon, lat = weight$lat, + filled.continents = FALSE, title_scale = 0.8, + bar_label_scale = 1.3, axes_label_scale = 1.2, + intylat = 2, intxlon = 2, degree_sym = TRUE, + toptitle = 'December Weights nf 4') +dev.off() +agg_png(paste0(dir_output, "Slope.png"), + width = 700, height = 700, units = 'px',res = 144) +plot(1:12, slope_plot, type = 'b', main = 'Slope', pch = 16, xlab = 'month', + ylab = 'Slope', bty = 'n', cex.main = 1.5, cex.lab = 1.3, cex = 1.5) +lines(12, slope_plot[12], type = 'p', col = 'red', pch = 16, cex = 1.3) +dev.off() + +# Plot ForecastPDF +library(abind) +#obsts <- MeanDims(obs$data, c('lat', 'lon'), na.rm = T) +#print(quantile(obsts, c(0.1, 0.3, 0.6, 0.9), na.rm = T)) +#expts <- MeanDims(exp$data, c('lat', 'lon'), na.rm = T) +#exp.qmts <- MeanDims(exp.qm$data, c('lat', 'lon'), na.rm = T) +print("Quantiles gridpoint") +print(quantile(as.vector(exp.qm$data[1,,,11,5,4,2]), c(0.1,0.3,0.6,0.9))) +data <- rbind(exp$data[1, ,1,11,5,4,2], exp.qm$data[1,,1,11,5,4,2]) +empty <- array(NA, c(dataset = 2, members = 225)) +names(dim(data)) <- names(dim(empty)) +data <- abind(data, empty, along = 2) +names(dim(data)) <- names(dim(empty)) +data <- abind(data, fs$data[1,,1,11,15,18 ,2], along = 1) +names(dim(data)) <- names(dim(empty)) +print(dim(data)) +print(quantile(obs$data[1,,,11,5,4,2])) +agg_png(paste0(dir_output, "/FiguresPDF_RF4_11December_GridPoint.png"), + width = 1000, height = 1000, units = 'px', res = 144) +PlotForecastPDF(data, tercile.limits = c(0.02, 1.17), + extreme.limits = c(0.0155555, 7.75), color.set = 'hydro', + var.name = "Precipitation (mm)", add.ensmemb = 'no', + title = "Forecasts issued on Nov 1993 for 11th December 1993", + fcst.names = c("SEAS5", "Bias Adjusted", "Downscaled nf 4")) +dev.off() + +obsts <- MeanDims(obs$data, c('lat', 'lon'), na.rm = T) +print(quantile(obsts, c(0.1, 0.3, 0.6, 0.9), na.rm = T)) +expts <- MeanDims(exp$data, c('lat', 'lon'), na.rm = T) +exp.qmts <- MeanDims(exp.qm$data, c('lat', 'lon'), na.rm = T) +empty <- array(NA, c(dataset = 2, member = 225, sdate = 26, ftime = 31, monthly = 8)) +data <- abind(expts, exp.qmts, along = 1) +names(dim(data)) <- names(dim(expts)) +data <- abind(data, empty, along = 2) +names(dim(data)) <- names(dim(expts)) +fsts <- MeanDims(fs$data, c('lat', 'lon'), na.rm = T) +data <- abind(data, fsts, along = 1) +names(dim(data)) <- c('dataset', 'members', 'sdate', 'ftime', 'monthly') +agg_png(paste0(dir_output, "/FiguresPDF_RF4_11December_Areave.png"), + width = 1400, height = 800, units = 'px', res = 144) +PlotForecastPDF(data[,,1,11,2], tercile.limits = c(0.67, 2.5), + extreme.limits = c(0.09, 7.3), color.set = 'hydro', + var.name = "Precipitation (mm)", + title = "Forecasts issued on Nov 1993 for 11th December 1993", + fcst.names = c("SEAS5", "Bias Adjusted", "Downscaled nf 4")) +dev.off() + + + + + + diff --git a/inst/doc/UseCase3_data_preparation_SCHEME_model.R b/inst/doc/UseCase3_data_preparation_SCHEME_model.R new file mode 100644 index 0000000000000000000000000000000000000000..ada24ef2c762ca78eef563a2fca2c97a75a44903 --- /dev/null +++ b/inst/doc/UseCase3_data_preparation_SCHEME_model.R @@ -0,0 +1,512 @@ +# Author: Bert Van Schaeybroeck +# Use Case 3: Seasonal forecasts for a river flow +# ----------------------------------------------- +rm(list = ls()) +library(CSTools) +library(s2dverification) +library(CSTools) +library(ClimProjDiags) + +#SETUP PARAMETERS (TO FIX BEFORE RUNNING SCRIPT): +#------------------------------------------------ +var.to.use <- "prlr" #which variable to correct (prlr=precipitation, tasmin, tasmax, tas) + +init.yr <- 1993 #initial year (for ECMWF Sys5 = 1993) +end.yr <- 2019 #end year (for ECMWF Sys5 = 2019) +amt.ftime <- 214 +n.cores.to.use <- 20 +use.chirps <- T +eval.method.to.use <- "leave-one-out" +domain.high.res <- "greece_high_res" +domain.low.res <- "greece_low_res" +sdate.mon.to.use <- 5 #Month of startdate for ECMWF Sys5 possibilities are 5 (May) or 11 (November) +sdate.day.to.use <- 1 #Day of startdate for ECMWF Sys5 only possibility is 1 +make.plot.msk <- F #mask to indicate if figures need to be made based on the output of the Analog function. + +#LOCAL PARAMETERS (to be adjusted for each working system) +#--------------------------------------------------------- +dir.rdata <- "/mnt/HDS_URCLIM/URCLIM/bertvs/medscope/data/greece_rdata/" + +#dir.scripts <- "/mnt/netapp/home/bertvs/ARCHIVE_bertvs/R/medscope/D3.1/" +dir.scratch <- "/scratch/bertvs/" +#Base path for C3S forecast (experiment) dataset: +dir.c3s <- "/mnt/HDS_URCLIM/URCLIM/bertvs/medscope/data/" +#Base path for ERA5 reference or obs dataset: +dir.era5 <- "/mnt/HDS_BREGILABEPOC/BREGILABEPOC/era5/europe/daily/per_mon/" +#Base path for CHIRPS (reference or obs) rainfall dataset: +dir.chirps <- "/mnt/HDS_MEDYCLIM/MEDYCLIM/PREDANTAR/climate_data/obs/chirps/" +dir.chirps.low.res <- paste0(dir.chirps, "/", domain.low.res, "/per_mon/") +dir.chirps.high.res <- paste0(dir.chirps, "/", domain.high.res, "/per_mon/") + + +#AUXILIARY FUNCTIONS +#------------------- +set.msk <- function(x, msk, const){ + x[msk] = const + return(x) +} + +#FIXED PARAMETERS: +#----------------- +greece.coor.vec <- c( + lonmin = 18.975, + lonmax = 24.025, + latmin = 37.975, + latmax = 43.025) +greece.coor.lst <- list( + lon.min = 18.975, + lon.max = 24.025, + lat.min = 37.975, + lat.max = 43.025) +coor.to.use <- greece.coor.lst + +europe.coor <- list( + lon.min = -27, + lon.max = 45, + lat.min = 33, + lat.max = 73.5) + +#Large-scale pressure field metadata (necessary for analogs) +var.msl <- "mean_sea_level_pressure" +nc.var.name.msl <- "msl" + +#Depending on the variable loaded, different datasets and metadata are used +if(var.to.use == "prlr"){ #Precipitation + var.era5 <- "total_precipitation" + time.era5 <- "daily" + nc.var.name.era5 <- "tp" + var.chirps <- "precip" + time.chirps <- "daily" + nc.var.name.chirps <- "precip" + cal.meth.to.use <- "bias" #method for bias calibration + chirps.low.res.daily <- list( + name = "chirps_low_res", + path = paste0(dir.chirps.low.res, "chirps-v2.0.$YEAR$.days_greece_low_res-$MONTH$.nc"), + nc_var_name = nc.var.name.chirps) + chirps.high.res.daily <- list( + name = "chirps_high_res", + path = paste0(dir.chirps.high.res, + "chirps-v2.0.$YEAR$.days_greece_high_res-$MONTH$.nc"), + nc_var_name = nc.var.name.chirps) + #unit conversions + mul.cor.era5 <- 1000 * 24 + add.cor.era5 <- 0 + mul.cor.chirps <- 1 + add.cor.chirps <- 0 + mul.cor.exp <- 1000 * 3600 * 24 + add.cor.exp <- 0 + if(use.chirps){ + add.cor.obs <- add.cor.chirps + mul.cor.obs <- mul.cor.chirps + } else { + add.cor.obs <- add.cor.era5 + mul.cor.obs <- mul.cor.er5 + } +} else if(var.to.use == "tas"){ + var.era5 <- "2m_temperature" + time.era5 <- "daily" + nc.var.name.era5 <- "t2m" + cal.meth.to.use <- "mse_min" + + #unit conversions + mul.cor.era5 <- 0 + add.cor.era5 <- 0 + mul.cor.exp <- 0 + add.cor.exp <- 0 +} else if(var.to.use == "tasmin"){ + var.era5 <- "2m_temperature" + time.era5 <- "daily_min" + nc_var_name.era5 <- "t2m" + cal.meth.to.use <- "mse_min" #method for bias calibration + + #unit conversions + mul.cor.era5 <- 0 + add.cor.era5 <- 0 + mul.cor.exp <- 0 + add.cor.exp <- 0 +} else if(var.to.use == "tasmax"){ + var.era5 <- "2m_temperature" + time.era5 <- "daily_max" + nc_var_name.era5 <- "t2m" + cal.meth.to.use <- "mse_min" #method for bias calibration + + #unit conversions + mul.cor.era5 <- 0 + add.cor.era5 <- 0 + mul.cor.exp <- 0 + add.cor.exp <- 0 +} + + +#Experiment path specification: +ecmwf.s5.daily <- list( + name = "ecmwfS5", + path = paste0(dir.c3s, + "C3S/$EXP_NAME$/$STORE_FREQ$/$VAR_NAME$/", + "$VAR_NAME$_$START_DATE$.nc")) +#Reference or obs path specifications (ERA5 data available over Europe): +era5.daily <- list(name = "era5", + path = paste0(dir.era5, "era5-", time.era5, + "-europe-", var.era5, "-$YEAR$-$MONTH$.nc"), + nc_var_name = nc.var.name.era5) +#Reference or obs path specifications for pressure field (ERA5 data available over Europe): +msl.era5.daily <- list(name = "msl", + path = paste0(dir.era5, "era5-", time.era5, + "-europe-", var.msl, "-$YEAR$-$MONTH$.nc"), + nc_var_name = nc.var.name.msl) + +#UNIVERSAL PARAMETERS: +#--------------------- +amt.mon.per.yr <- 12 +amt.day.per.mon <- 31 +sdate.day.str <- formatC(sdate.day.to.use, width = 2, flag = "0") +sdate.mon.str <- formatC(sdate.mon.to.use, width = 2, flag = "0") +day.lst <- formatC(seq(1, amt.day.per.mon), width = 2, flag = "0") +yr.lst <- seq(init.yr, end.yr) +amt.yr <- length(yr.lst) +sdate.lst <- paste0(yr.lst, sdate.mon.str, sdate.day.str ) + + +#START +#----- + +#1. LOAD THE DATA +#---------------- + +if(use.chirps){ + obs.set.to.use <- chirps.low.res.daily +} else { + obs.set.to.use <- era5.daily +} + +#Load mean sea level pressure field from ERA5 (no need to set the units) + +file.to.load <- paste0(dir.rdata, "msl_all.RData") +if(file.exists(file.to.load)){ + load(file.to.load, verbose = T) +} else { + msl.all <- CST_Load( + var = "psl", #nc.var.name.msl, + obs = list(msl.era5.daily), + exp = list(ecmwf.s5.daily), + nmember = NULL, + sdates = sdate.lst, + lonmin = europe.coor$lon.min, + lonmax = europe.coor$lon.max, + latmin = europe.coor$lat.min, + latmax = europe.coor$lat.max, + output = "lonlat", + storefreq = "daily", + nprocs = n.cores.to.use) + save(file = file.to.load, msl.all) +} + +#Data manipulation: first split lead times per month. Then merge all data per month and all sdates. +#This merged dataset will be used to calibrate (per month) and find analogs (per month). +obs.msl.eur.split <- CST_SplitDim(msl.all$obs, split_dim = c("ftime")) +exp.msl.eur.split <- CST_SplitDim(msl.all$exp, split_dim = c("ftime")) +obs.msl.eur.merge <- CST_MergeDims( + obs.msl.eur.split, + merge_dims = c("sdate", "ftime"), + rename_dim = "sdate") +exp.msl.eur.merge <- CST_MergeDims( + exp.msl.eur.split, + merge_dims = c("sdate", "ftime"), + rename_dim = "sdate") + +obs.msl.eur.merge.an <- CST_Anomaly(exp = obs.msl.eur.merge, dim_anom = 3) +exp.msl.eur.merge.an <- CST_Anomaly(exp = exp.msl.eur.merge, dim_anom = 3) + + +#Load observational and forecast set of variable that needs to be calibrated and downscaled: +file.to.load <- paste0(dir.rdata, "data_all.RData") +if(file.exists(file.to.load)){ + load(file.to.load, verbose = T) +} else { + data.all <- CST_Load( + var = var.to.use, + obs = list(obs.set.to.use), + exp = list(ecmwf.s5.daily), + nmember = NULL, + sdates = sdate.lst, + lonmin = coor.to.use$lon.min, + lonmax = coor.to.use$lon.max, + latmin = coor.to.use$lat.min, + latmax = coor.to.use$lat.max, + output = "lonlat", + storefreq = "daily", + nprocs = n.cores.to.use) + save(file = file.to.load, data.all) +} +#Set the units: +data.all$obs$data <- data.all$obs$data * mul.cor.obs + add.cor.obs +data.all$exp$data <- data.all$exp$data * mul.cor.exp + add.cor.exp + +#Data manipulation: first split lead times per month. Then merge all data per month and all sdates. +#This merged dataset will be used to calibrate (per month) and find analogs (per month). +obs.split <- CST_SplitDim(data.all$obs, split_dim = c("ftime")) +exp.split <- CST_SplitDim(data.all$exp, split_dim = c("ftime")) +obs.merge <- CST_MergeDims( + obs.split, + merge_dims = c("sdate", "ftime"), + rename_dim = "sdate") +exp.merge <- CST_MergeDims( + exp.split, + merge_dims = c("sdate", "ftime"), + rename_dim = "sdate") + +#Calibrate the exp data (per month) +cal.merge <- CST_Calibration( + exp = exp.merge, + obs = obs.merge, + cal.method = cal.meth.to.use, + eval.method = eval.method.to.use) +cal.merge$data[cal.merge$data < 0] <- 0 + +#LOAD HIGH RES CHIRPS DATA +file.to.load <- paste0(dir.rdata, "obs_high_res.RData") +if(file.exists(file.to.load)){ + load(file.to.load, verbose = T) +} else { + obs.high.res <- CST_Load(var = var.to.use, + obs = list(chirps.high.res.daily), + exp = NULL, + sdates = sdate.lst, + nmember = 1, + leadtimemax = amt.ftime, + sampleperiod = 1, + lonmin = coor.to.use$lon.min, + lonmax = coor.to.use$lon.max, + latmin = coor.to.use$lat.min, + latmax = coor.to.use$lat.max, + output = "lonlat", + storefreq = "daily", + nprocs = n.cores.to.use) + save(file = file.to.load, obs.high.res) +} +#set the units +obs.high.res$data <- obs.high.res$data * mul.cor.chirps + + add.cor.chirps +#split per month +obs.high.res.split <- CST_SplitDim( + obs.high.res, + split_dim = c("ftime")) +#merge lead times and sdates +obs.high.res.merge <- CST_MergeDims( + obs.high.res.split, + merge_dims = c("sdate", "ftime"), + rename_dim = "sdate") + +#LOAD LOW RES CHIRPS DATA +file.to.load <- paste0(dir.rdata, "obs_low_res.RData") +if(file.exists(file.to.load)){ + load(file.to.load, verbose = T) +} else { + obs.low.res <- CST_Load(var = var.to.use, + obs = list(chirps.low.res.daily), + exp = NULL, + sdates = sdate.lst, + nmember = 1, + leadtimemax = amt.ftime, + sampleperiod = 1, + lonmin = coor.to.use$lon.min, + lonmax = coor.to.use$lon.max, + latmin = coor.to.use$lat.min, + latmax = coor.to.use$lat.max, + output = "lonlat", + storefreq = "daily", + nprocs = n.cores.to.use) + save(file = file.to.load, obs.low.res) +} +#set units +obs.low.res$data <- obs.low.res$data * mul.cor.chirps + + add.cor.chirps +#split per month +obs.low.res.split <- CST_SplitDim( + obs.low.res, + split_dim = c("ftime")) +#merge lead times and sdates +obs.low.res.merge <- CST_MergeDims( + obs.low.res.split, + merge_dims = c("sdate", "ftime"), + rename_dim = "sdate") + + +#2. PROCESS THE DATA +#------------------- + +#amount of ensemble members from experiment. For ECMWF Sys5 it is 25: +amt.mbr <- as.numeric(dim(cal.merge$data)["member"]) +lon.low.res <- as.vector(cal.merge$lon) +lat.low.res <- as.vector(cal.merge$lat) +lon.high.res <- as.vector(obs.high.res$lon) +lat.high.res <- as.vector(obs.high.res$lat) +lon.eur <- as.vector(obs.msl.eur.merge.an$lon) +lat.eur <- as.vector(obs.msl.eur.merge.an$lat) + +#amount of lead times in months. For ECMWF Sys5 it is 7: +amt.lead.mon <- as.numeric(dim(cal.merge$data)["monthly"]) +mon.seq.tmp <- seq(sdate.mon.to.use, sdate.mon.to.use + amt.lead.mon - 1) +mon.seq.tmp <- ((mon.seq.tmp - 1) %% amt.mon.per.yr) + 1 +lead.mon.lst <- formatC(mon.seq.tmp, width = 2, flag = "0") +#amount of starting days from experiment. For ECMWF Sys5 it is 837: +amt.sdate <- as.numeric(dim(cal.merge$data)["sdate"]) + +sub.time <- outer( + as.vector(t(outer(yr.lst, day.lst, paste, sep="-"))), + lead.mon.lst, + paste, sep = "-") +#This step is necessary to set the non-existent dates to NA +sub.time <- format(as.Date(sub.time, format("%Y-%d-%m")), "%Y-%m-%d") +dim(sub.time) <- c(sdate = amt.yr * amt.day.per.mon, time = amt.lead.mon) + +cal.high.res.merge <- obs.high.res.merge +cal.high.res.merge$data[] <- NA + +#Determine spatial points with all obs.high.res.merge (CHIRPS) data equal to NA. These are the points over sea. +is.na.high.res.obs <- apply( + obs.high.res.merge$data, + MARGIN = c(4, 5), + FUN = function(x){all(is.na(x))}) +#Determine spatial points with all obs.low.res.merge (CHIRPS) data equal to NA. These are the points over sea. +is.na.low.res.obs <- apply( + obs.low.res.merge$data, + MARGIN = c(4, 5), + FUN = function(x){all(is.na(x))}) + +#Set all calibrated exp data (cal.merge) equal to NA at the sea point. +cal.merge.tmp = Apply( + data = list(x = cal.merge$data), + target_dims = list(x = c("lat", "lon")), + fun = set.msk, + msk = is.na.low.res.obs, + const = 0, + output_dims = list(c("lat", "lon")) + )$output1 +dex.match <- match(names(dim(cal.merge$data)), names(dim(cal.merge.tmp))) +cal.merge$data <- aperm(cal.merge.tmp, dex.match) +rm(cal.merge.tmp) + +#2. PROCESS THE DATA +#------------------- + +i.dataset <- 1 +i.mbr.obs <- 1 +for(i.mbr in seq(1, amt.mbr)){ + for(i.mon in seq(1, amt.lead.mon)){ + for(i.sdate in seq(1, amt.sdate)){ + #i.mbr <- 1 + #i.mon = 1 + #i.sdate = 24 + date.to.use <- sub.time[ i.sdate, i.mon] + date.an.lst <- sub.time[ , i.mon] + cat("i.mbr = ", i.mbr, ", i.mon =", i.mon, ", i.sdate = ", + i.sdate, "date: ", date.to.use,"\n") + + + #Extract the (calibrated) forecast that you want to downscale: + exp.low.res.tmp <- exp.merge$data[i.dataset, i.mbr, i.sdate, , , i.mon] + cal.low.res.tmp <- cal.merge$data[i.dataset, i.mbr, i.sdate, , , i.mon] + #Extract the large-scale pressure field of that day + exp.msl.eur.tmp <- exp.msl.eur.merge.an$data[i.dataset, i.mbr, i.sdate, , , i.mon] + + #Extract all observations that will be used to find analogs + obs.msl.eur.tmp <- obs.msl.eur.merge.an$data[i.dataset, i.mbr.obs, , , , i.mon]#-i.sdate + obs.low.res.tmp <- obs.low.res.merge$data[i.dataset, i.mbr.obs, , , , i.mon] #-i.sdate + obs.high.res.tmp <- obs.high.res.merge$data[i.dataset, i.mbr.obs, , , , i.mon] #-i.sdate + names(dim(obs.high.res.tmp)) <- c("time", "lat", "lon") + names(dim(obs.low.res.tmp)) <- c("time", "lat", "lon") + names(dim(obs.msl.eur.tmp)) <- c("time", "lat", "lon") + if(!is.na(date.to.use) & !all(is.na(cal.low.res.tmp))){ + obs.low.res.tmp[is.na(obs.low.res.tmp)] <- 0 + + res <- Analogs( + expL = exp.msl.eur.tmp, + obsL = obs.msl.eur.tmp, + time_obsL = date.an.lst, + obsVar = obs.low.res.tmp, + expVar = exp.low.res.tmp, + lonVar = lon.low.res, + latVar = lat.low.res, + lonL = lon.eur, + latL = lat.eur, + region = greece.coor.vec, + criteria = "Local_dist", + time_expL = date.to.use, + excludeTime = date.to.use, + AnalogsInfo = T, + nAnalogs = 1000) + + + if(make.plot.msk){ + corr.date <- as.character(res$dates[1]) #select the date of the most + corr.dex <- which(date.an.lst == corr.date) + + #The following figure shows the uncalibrated raw model field (analogous to Fig. 9a) + file.fig <- paste0("mbr_", i.mbr, "_mon_", i.mon, + "_sdate_", date.to.use, "_exp.low.res.pdf") + pdf(file = file.fig) + PlotEquiMap( + exp.low.res.tmp[ , ], + lon = obs.low.res.merge$lon, + lat = obs.low.res.merge$lat, + filled.continents = F, + intylat = 2, + intxlon = 2, + title_scale = 0.7, #bar_limits = c(0, 60), + units = "precipitation (mm)") + dev.off() + + #The following figure includes the calibrated model field (analogous to Fig. 9b) + file.fig <- paste0("mbr_", i.mbr, "_mon_", i.mon, + "_sdate_", date.to.use, "_cal.low.res.pdf") + pdf(file = file.fig) + PlotEquiMap( + cal.low.res.tmp, + lon = obs.low.res.merge$lon, + lat = obs.low.res.merge$lat, + filled.continents = F, + intylat = 2, + intxlon = 2, + title_scale = 0.7, #bar_limits = c(0, 60), + units = "precipitation (mm)") + + #The following figure includes the analog upscaled field (analogous to Fig. 9c) + file.fig <- paste0("mbr_", i.mbr, "_mon_", i.mon, + "_sdate_", date.to.use, "_obs.low.res.pdf") + pdf(file = file.fig) + PlotEquiMap( + obs.low.res.tmp[corr.dex, , ], + lon = obs.low.res.merge$lon, + lat = obs.low.res.merge$lat, + filled.continents = F, + intylat = 2, + intxlon = 2, + title_scale = 0.7, #bar_limits = c(0, 60), + units = "precipitation (mm)") + dev.off() + + #The following figure includes the analog field (analogous to Fig. 9d) + file.fig <- paste0("mbr_", i.mbr, "_mon_", i.mon, + "_sdate_", date.to.use, "_obs.high.res.pdf") + pdf(file = file.fig) + PlotEquiMap( + obs.high.res.tmp[corr.dex, , ], + lon = obs.high.res.merge$lon, + lat = obs.high.res.merge$lat, + filled.continents = F, + intylat = 2, + intxlon = 2, + title_scale = 0.7, #bar_limits = c(0, 60), + units = "precipitation (mm)") + dev.off() + + + } + } + } + } +} + + + diff --git a/inst/doc/launch_UseCase2_PrecipitationDownscaling_RF4.sh b/inst/doc/launch_UseCase2_PrecipitationDownscaling_RF4.sh new file mode 100644 index 0000000000000000000000000000000000000000..2e1dd93945d0794a4451584920c45e8791d46c38 --- /dev/null +++ b/inst/doc/launch_UseCase2_PrecipitationDownscaling_RF4.sh @@ -0,0 +1,10 @@ +#!/bin/bash +#BSUB -W 6:00 +#BSUB -n 16 +#BSUB -M 7000 +#BSUB -J RainFARM_Downsc +#BSUB -o /esarchive/scratch/nperez/CSTools_manuscript/v20210603/lsf-%J.out +#BSUB -e /esarchive/scratch/nperez/CSTools_manuscript/v20210603/lsf-%J.err + +module load R +Rscript /esarchive/scratch/nperez/git/cstools/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R diff --git a/man/Analogs.Rd b/man/Analogs.Rd index 746ebdd1d8e48fac7a7b21b68a1f029626ec71b8..fc26a5523fe0dd78c4755b03f106118d99b193b3 100644 --- a/man/Analogs.Rd +++ b/man/Analogs.Rd @@ -9,6 +9,8 @@ Analogs( obsL, time_obsL, time_expL = NULL, + lonL = NULL, + latL = NULL, expVar = NULL, obsVar = NULL, criteria = "Large_dist", @@ -41,6 +43,10 @@ in the format "dd/mm/yyyy". Reference time to search for analogs.} dimensions in expL) of character string(s) indicating the date(s) of the experiment in the format "dd/mm/yyyy". Time(s) to find the analogs.} +\item{lonL}{a vector containing the longitude of parameter 'expL'.} + +\item{latL}{a vector containing the latitude of parameter 'expL'.} + \item{expVar}{an array of N named dimensions containing the experimental field on the local scale, usually a different variable to the parameter 'expL'. If it is not NULL (by default, NULL), the returned field by this @@ -167,21 +173,21 @@ downscale_field <- Analogs(expL = expSLP, obsL = obsSLP, obsVar = obs.pr, region=c(lonmin = -1 ,lonmax = 2, latmin = 30, latmax = 33) Local_scale <- Analogs(expL = expSLP, obsL = obsSLP, time_obsL = time_obsSLP, obsVar = obs.pr, criteria = "Local_dist", - lonVar = seq(-1, 5, 1.5),latVar = seq(30, 35, 1.5), + lonL = seq(-1, 5, 1.5),latL = seq(30, 35, 1.5), region = region,time_expL = "01-10-2000", nAnalogs = 10, AnalogsInfo = TRUE) # Example 6: list of best analogs using criteria 'Local_dist' and 2 Local_scale <- Analogs(expL = expSLP, obsL = obsSLP, time_obsL = time_obsSLP, - criteria = "Local_dist", lonVar = seq(-1, 5, 1.5), - latVar = seq(30, 35, 1.5), region = region, + criteria = "Local_dist", lonL = seq(-1, 5, 1.5), + latL = seq(30, 35, 1.5), region = region, time_expL = "01-10-2000", nAnalogs = 5, AnalogsInfo = TRUE) # Example 7: Downscaling using Local_dist criteria Local_scale <- Analogs(expL = expSLP, obsL = obsSLP, time_obsL = time_obsSLP, - criteria = "Local_dist", lonVar = seq(-1, 5, 1.5), - latVar = seq(30, 35, 1.5), region = region, + criteria = "Local_dist", lonL = seq(-1, 5, 1.5), + latL = seq(30, 35, 1.5), region = region, time_expL = "01-10-2000", nAnalogs = 10, AnalogsInfo = FALSE) @@ -190,14 +196,16 @@ exp.pr <- c(rnorm(1:20) * 0.001) dim(exp.pr) <- dim(expSLP) Local_scalecor <- Analogs(expL = expSLP, obsL = obsSLP,time_obsL = time_obsSLP, obsVar = obs.pr, expVar = exp.pr, - criteria = "Local_cor", lonVar = seq(-1, 5, 1.5), - time_expL = "01-10-2000", latVar = seq(30, 35, 1.5), + criteria = "Local_cor", lonL = seq(-1, 5, 1.5), + time_expL = "01-10-2000", latL = seq(30, 35, 1.5), + lonVar = seq(-1, 5, 1.5), latVar = seq(30, 35, 1.5), nAnalogs = 8, region = region, AnalogsInfo = FALSE) # same but without imposing nAnalogs,so nAnalogs will be set by default as 10 Local_scalecor <- Analogs(expL = expSLP, obsL = obsSLP,time_obsL = time_obsSLP, obsVar = obs.pr, expVar = exp.pr, - criteria = "Local_cor", lonVar = seq(-1,5,1.5), - time_expL = "01-10-2000", latVar=seq(30, 35, 1.5), + lonVar = seq(-1, 5, 1.5), latVar = seq(30, 35, 1.5), + criteria = "Local_cor", lonL = seq(-1,5,1.5), + time_expL = "01-10-2000", latL =seq(30, 35, 1.5), region = region, AnalogsInfo = TRUE) #'Example 9: List of best analogs in the three criterias Large_dist, @@ -206,11 +214,12 @@ Large_scale <- Analogs(expL = expSLP, obsL = obsSLP, time_obsL = time_obsSLP, nAnalogs = 7, AnalogsInfo = TRUE) Local_scale <- Analogs(expL = expSLP, obsL = obsSLP, time_obsL = time_obsSLP, time_expL = "01-10-2000", criteria = "Local_dist", - lonVar = seq(-1, 5, 1.5), latVar = seq(30, 35, 1.5), + lonL = seq(-1, 5, 1.5), latL = seq(30, 35, 1.5), nAnalogs = 7,region = region, AnalogsInfo = TRUE) Local_scalecor <- Analogs(expL = expSLP, obsL = obsSLP,time_obsL = time_obsSLP, obsVar = obsSLP, expVar = expSLP, time_expL = "01-10-2000",criteria = "Local_cor", + lonL = seq(-1, 5, 1.5), latL = seq(30, 35, 1.5), lonVar = seq(-1, 5, 1.5), latVar = seq(30, 35, 1.5), nAnalogs = 7,region = region, AnalogsInfo = TRUE) diff --git a/man/BiasCorrection.Rd b/man/BiasCorrection.Rd new file mode 100644 index 0000000000000000000000000000000000000000..8b60a6de655fd4261a51bf72ad02f99ac2756498 --- /dev/null +++ b/man/BiasCorrection.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_BiasCorrection.R +\name{BiasCorrection} +\alias{BiasCorrection} +\title{Bias Correction based on the mean and standard deviation adjustment} +\usage{ +BiasCorrection(exp, obs, exp_cor = NULL, na.rm = FALSE) +} +\arguments{ +\item{exp}{a multidimensional array with named dimensions containing the seasonal forecast experiment data with at least 'member' and 'sdate' dimensions.} + +\item{obs}{a multidimensional array with named dimensions containing the observed data with at least 'sdate' dimension.} + +\item{exp_cor}{a multidimensional array with named dimensions containing the seasonl forecast experiment to be corrected. If it is NULL, the 'exp' forecast will be corrected.} + +\item{na.rm}{a logical value indicating whether missing values should be stripped before the computation proceeds, by default it is set to FALSE.} +} +\value{ +an object of class \code{s2dv_cube} containing the bias corrected forecasts in the element called \code{$data} with the same dimensions of the experimental data. +} +\description{ +This function applies the simple bias adjustment technique described in Torralba et al. (2017). The adjusted forecasts have an equivalent standard deviation and mean to that of the reference dataset. +} +\examples{ + +# Example +# Creation of sample s2dverification objects. These are not complete +# s2dverification objects though. The Load function returns complete objects. +mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7) +dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7) +obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7) +dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7) +a <- BiasCorrection(exp = mod1, obs = obs1) +str(a) +} +\references{ +Torralba, V., F.J. Doblas-Reyes, D. MacLeod, I. Christel and M. Davis (2017). Seasonal climate prediction: a new source of information for the management of wind energy resources. Journal of Applied Meteorology and Climatology, 56, 1231-1247, doi:10.1175/JAMC-D-16-0204.1. (CLIM4ENERGY, EUPORIAS, NEWA, RESILIENCE, SPECS) +} +\author{ +Verónica Torralba, \email{veronica.torralba@bsc.es} +} diff --git a/man/CST_BiasCorrection.Rd b/man/CST_BiasCorrection.Rd index 55c325a2db6dc53f8b315c4c0afc718746f79b1a..adfc2798ac9c9ddfa3cb942a0c2faf69233bfaed 100644 --- a/man/CST_BiasCorrection.Rd +++ b/man/CST_BiasCorrection.Rd @@ -4,13 +4,15 @@ \alias{CST_BiasCorrection} \title{Bias Correction based on the mean and standard deviation adjustment} \usage{ -CST_BiasCorrection(exp, obs, na.rm = FALSE) +CST_BiasCorrection(exp, obs, exp_cor = NULL, na.rm = FALSE) } \arguments{ \item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}} \item{obs}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed data in the element named \code{$data}.} +\item{exp_cor}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonl forecast experiment to be corrected. If it is NULL, the 'exp' forecast will be corrected.} + \item{na.rm}{a logical value indicating whether missing values should be stripped before the computation proceeds, by default it is set to FALSE.} } \value{ diff --git a/man/CST_DynBiasCorrection.Rd b/man/CST_DynBiasCorrection.Rd new file mode 100644 index 0000000000000000000000000000000000000000..e2467e72bbc98eee7b45c0b5bd6192506bf22fcf --- /dev/null +++ b/man/CST_DynBiasCorrection.Rd @@ -0,0 +1,98 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_DynBiasCorrection.R +\name{CST_DynBiasCorrection} +\alias{CST_DynBiasCorrection} +\title{Performing a Bias Correction conditioned by the dynamical +properties of the data.} +\usage{ +CST_DynBiasCorrection( + exp, + obs, + method = "QUANT", + wetday = FALSE, + proxy = "dim", + quanti, + ncores = NULL +) +} +\arguments{ +\item{exp}{an s2v_cube object with the experiment data} + +\item{obs}{an s2dv_cube object with the reference data} + +\item{method}{a character string indicating the method to apply bias +correction among these ones: "PTF","RQUANT","QUANT","SSPLIN"} + +\item{wetday}{logical indicating whether to perform wet day correction +or not OR a numeric threshold below which all values are set to zero (by +default is set to 'FALSE').} + +\item{proxy}{a character string indicating the proxy for local dimension +'dim' or inverse of persistence 'theta' to apply the dynamical +conditioned bias correction method.} + +\item{quanti}{a number lower than 1 indicating the quantile to perform +the computation of local dimension and theta} + +\item{ncores}{The number of cores to use in parallel computation} +} +\value{ +dynbias an s2dvcube object with a bias correction performed +conditioned by local dimension 'dim' or inverse of persistence 'theta' +} +\description{ +This function perform a bias correction conditioned by the +dynamical properties of the dataset. This function internally uses the functions +'Predictability' to divide in terciles the two dynamical proxies +computed with 'CST_ProxiesAttractor'. A bias correction +between the model and the observations is performed using the division into +terciles of the local dimension 'dim' and inverse of the persistence 'theta'. +For instance, model values with lower 'dim' will be corrected with observed +values with lower 'dim', and the same for theta. The function gives two options +of bias correction: one for 'dim' and/or one for 'theta' +} +\examples{ +# example 1: simple data s2dvcube style +set.seed(1) +expL <- rnorm(1:2000) +dim (expL) <- c(time =100,lat = 4, lon = 5) +obsL <- c(rnorm(1:1980),expL[1,,]*1.2) +dim (obsL) <- c(time = 100,lat = 4, lon = 5) +time_obsL <- paste(rep("01", 100), rep("01", 100), 1920 : 2019, sep = "-") +time_expL <- paste(rep("01", 100), rep("01", 100), 1929 : 2019, sep = "-") +lon <- seq(-1,5,1.5) +lat <- seq(30,35,1.5) +# qm=0.98 # too high for this short dataset, it is possible that doesn't +# get the requirement, in that case it would be necessary select a lower qm +# for instance qm=0.60 +expL <- s2dv_cube(data = expL, lat = lat, lon = lon, + Dates = list(start = time_expL, end = time_expL)) +obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, + Dates = list(start = time_obsL, end = time_obsL)) +# to use DynBiasCorrection +dynbias1 <- DynBiasCorrection(exp = expL$data, obs = obsL$data, proxy= "dim", + quanti = 0.6) +# to use CST_DynBiasCorrection +dynbias2 <- CST_DynBiasCorrection(exp = expL, obs = obsL, proxy= "dim", + quanti = 0.6) + +} +\references{ +Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., +and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large +scale atmospheric predictability.Nature Communications, 10(1), 1316. +DOI = https://doi.org/10.1038/s41467-019-09305-8 " + +Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +Dynamical proxies of North Atlantic predictability and extremes. +Scientific Reports, 7-41278, 2017. +} +\author{ +Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} + +Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} + +Veronica Torralba, \email{veronica.torralba@cmcc.it} + +Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +} diff --git a/man/CST_ProxiesAttractor.Rd b/man/CST_ProxiesAttractor.Rd new file mode 100644 index 0000000000000000000000000000000000000000..ecf9b9da18bed146602da0ddb4f873cce452f9ec --- /dev/null +++ b/man/CST_ProxiesAttractor.Rd @@ -0,0 +1,51 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_ProxiesAttractor.R +\name{CST_ProxiesAttractor} +\alias{CST_ProxiesAttractor} +\title{Computing two dinamical proxies of the attractor in s2dv_cube.} +\usage{ +CST_ProxiesAttractor(data, quanti, ncores = NULL) +} +\arguments{ +\item{data}{a s2dv_cube object with the data to create the attractor. Must be a matrix with the timesteps in nrow +and the grids in ncol(dat(time,grids)} + +\item{quanti}{a number lower than 1 indicating the quantile to perform the computation of local dimension and theta} + +\item{ncores}{The number of cores to use in parallel computation} +} +\value{ +dim and theta +} +\description{ +This function computes two dinamical proxies of the attractor: +The local dimension (d) and the inverse of the persistence (theta) for an +'s2dv_cube' object. +These two parameters will be used as a condition for the computation of dynamical +scores to measure predictability and to compute bias correction conditioned by +the dynamics with the function DynBiasCorrection +Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in +} +\examples{ +# Example 1: Computing the attractor using simple s2dv data +attractor <- CST_ProxiesAttractor(data = lonlat_data$obs, quanti = 0.6) + +} +\references{ +Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019). +The hammam effect or how a warm ocean enhances large scale atmospheric predictability. +Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 " + +Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +Dynamical proxies of North Atlantic predictability and extremes. +Scientific Reports, 7-41278, 2017. +} +\author{ +Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} + +Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} + +Veronica Torralba, \email{veronica.torralba@cmcc.it} + +Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +} diff --git a/man/CST_SaveExp.Rd b/man/CST_SaveExp.Rd index ddd9164e1c9d88e284d4653d1b05941699e2adf4..92c972826018a806b67c84aa82644350a01130c6 100644 --- a/man/CST_SaveExp.Rd +++ b/man/CST_SaveExp.Rd @@ -5,7 +5,7 @@ \title{Save CSTools objects of class 's2dv_cube' containing experiments or observed data in NetCDF format} \usage{ -CST_SaveExp(data, destination = "./CST_Data") +CST_SaveExp(data, destination = "./CST_Data", extra_string = NULL) } \arguments{ \item{data}{an object of class \code{s2dv_cube}.} @@ -15,6 +15,8 @@ to save the data. NetCDF file for each starting date are saved into the folder tree: destination/experiment/variable/. By default the function creates and saves the data into the folder "CST_Data" in the working directory.} + +\item{extra_string}{a character string to be include as part of the file name, for instance, to identify member or realization. It would be added to the file name between underscore characters.} } \description{ This function allows to divide and save a object of class diff --git a/man/Calibration.Rd b/man/Calibration.Rd index f61a3cd3537844714533a899fa196125e9917319..7ac9cc2d11f272b7aa2bd05a07cca7c7c4843a42 100644 --- a/man/Calibration.Rd +++ b/man/Calibration.Rd @@ -48,7 +48,7 @@ Calibration( an array containing the calibrated forecasts with the same dimensions as the \code{exp} array. } \description{ -Four types of member-by-member bias correction can be performed. The \code{"bias"} method corrects the bias only, the \code{"evmos"} method applies a variance inflation technique to ensure the correction of the bias and the correspondence of variance between forecast and observation (Van Schaeybroeck and Vannitsem, 2011). The ensemble calibration methods \code{"mse_min"} and \code{"crps_min"} correct the bias, the overall forecast variance and the ensemble spread as described in Doblas-Reyes et al. (2005) and Van Schaeybroeck and Vannitsem (2015), respectively. While the \code{"mse_min"} method minimizes a constrained mean-squared error using three parameters, the \code{"crps_min"} method features four parameters and minimizes the Continuous Ranked Probability Score (CRPS). The \code{"rpc-based"} method adjusts the forecast variance ensuring that the ratio of predictable components (RPC) is equal to one, as in Eade et al. (2014). +Five types of member-by-member bias correction can be performed. The \code{"bias"} method corrects the bias only, the \code{"evmos"} method applies a variance inflation technique to ensure the correction of the bias and the correspondence of variance between forecast and observation (Van Schaeybroeck and Vannitsem, 2011). The ensemble calibration methods \code{"mse_min"} and \code{"crps_min"} correct the bias, the overall forecast variance and the ensemble spread as described in Doblas-Reyes et al. (2005) and Van Schaeybroeck and Vannitsem (2015), respectively. While the \code{"mse_min"} method minimizes a constrained mean-squared error using three parameters, the \code{"crps_min"} method features four parameters and minimizes the Continuous Ranked Probability Score (CRPS). The \code{"rpc-based"} method adjusts the forecast variance ensuring that the ratio of predictable components (RPC) is equal to one, as in Eade et al. (2014). Both in-sample or our out-of-sample (leave-one-out cross validation) calibration are possible. } diff --git a/man/DynBiasCorrection.Rd b/man/DynBiasCorrection.Rd new file mode 100644 index 0000000000000000000000000000000000000000..e6de373c814712165465b15d2f7802167b52982a --- /dev/null +++ b/man/DynBiasCorrection.Rd @@ -0,0 +1,83 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_DynBiasCorrection.R +\name{DynBiasCorrection} +\alias{DynBiasCorrection} +\title{Performing a Bias Correction conditioned by the dynamical +properties of the data.} +\usage{ +DynBiasCorrection( + exp, + obs, + method = "QUANT", + wetday = FALSE, + proxy = "dim", + quanti, + ncores = NULL +) +} +\arguments{ +\item{exp}{a multidimensional array with named dimensions with the +experiment data} + +\item{obs}{a multidimensional array with named dimensions with the +observation data} + +\item{method}{a character string indicating the method to apply bias +correction among these ones: +"PTF","RQUANT","QUANT","SSPLIN"} + +\item{wetday}{logical indicating whether to perform wet day correction +or not OR a numeric threshold below which all values are set to zero (by +default is set to 'FALSE').} + +\item{proxy}{a character string indicating the proxy for local dimension +'dim' or inverse of persistence 'theta' to apply the dynamical conditioned +bias correction method.} + +\item{quanti}{a number lower than 1 indicating the quantile to perform the +computation of local dimension and theta} + +\item{ncores}{The number of cores to use in parallel computation} +} +\value{ +a multidimensional array with named dimensions with a bias correction +performed conditioned by local dimension 'dim' or inverse of persistence 'theta' +} +\description{ +This function perform a bias correction conditioned by the +dynamical properties of the dataset. This function used the functions +'CST_Predictability' to divide in terciles the two dynamical proxies +computed with 'CST_ProxiesAttractor'. A bias correction +between the model and the observations is performed using the division into +terciles of the local dimension 'dim' and inverse of the persistence 'theta'. +For instance, model values with lower 'dim' will be corrected with observed +values with lower 'dim', and the same for theta. The function gives two options +of bias correction: one for 'dim' and/or one for 'theta' +} +\examples{ +expL <- rnorm(1:2000) +dim (expL) <- c(time =100,lat = 4, lon = 5) +obsL <- c(rnorm(1:1980),expL[1,,]*1.2) +dim (obsL) <- c(time = 100,lat = 4, lon = 5) +dynbias <- DynBiasCorrection(exp = expL, obs = obsL, method='QUANT', + proxy= "dim", quanti = 0.6) +} +\references{ +Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., +and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large +scale atmospheric predictability.Nature Communications, 10(1), 1316. +DOI = https://doi.org/10.1038/s41467-019-09305-8 " + +Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +Dynamical proxies of North Atlantic predictability and extremes. +Scientific Reports, 7-41278, 2017. +} +\author{ +Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} + +Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} + +Veronica Torralba, \email{veronica.torralba@cmcc.it} + +Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +} diff --git a/man/PlotCombinedMap.Rd b/man/PlotCombinedMap.Rd index c45d1afb594c362a6ecf8e4e641cf7223f8a0e1c..3d6661e14c45a3f9d5906125bce7d5a187e73d5a 100644 --- a/man/PlotCombinedMap.Rd +++ b/man/PlotCombinedMap.Rd @@ -16,8 +16,10 @@ PlotCombinedMap( col_unknown_map = "white", mask = NULL, col_mask = "grey", + dots = NULL, bar_titles = NULL, legend_scale = 1, + cex_bar_titles = 1.5, fileout = NULL, width = 8, height = 5, @@ -49,10 +51,19 @@ PlotCombinedMap( \item{col_mask}{Colour to be used for the superimposed mask (if specified in 'mask'). Takes the value 'grey' by default.} +\item{dots}{Array of same dimensions as 'var' or with dimensions +c(n, dim(var)), where n is the number of dot/symbol layers to add to the +plot. A value of TRUE at a grid cell will draw a dot/symbol on the +corresponding square of the plot. By default all layers provided in 'dots' +are plotted with dots, but a symbol can be specified for each of the +layers via the parameter 'dot_symbol'.} + \item{bar_titles}{Optional vector of character strings providing the titles to be shown on top of each of the colour bars.} \item{legend_scale}{Scale factor for the size of the colour bar labels. Takes 1 by default.} +\item{cex_bar_titles}{Scale factor for the sizes of the bar titles. Takes 1.5 by default.} + \item{fileout}{File where to save the plot. If not specified (default) a graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, bmp and tiff} \item{width}{File width, in the units specified in the parameter size_units (inches by default). Takes 8 by default.} diff --git a/man/Predictability.Rd b/man/Predictability.Rd new file mode 100644 index 0000000000000000000000000000000000000000..d37efcdcc4bf45c38ba5b4c175d919531e1e895f --- /dev/null +++ b/man/Predictability.Rd @@ -0,0 +1,76 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Predictability.R +\name{Predictability} +\alias{Predictability} +\title{Computing scores of predictability using two dynamical proxies +based on dynamical systems theory.} +\usage{ +Predictability(dim, theta, ncores = NULL) +} +\arguments{ +\item{dim}{An array of N named dimensions containing the local dimension as +the output of CST_ProxiesAttractor or ProxiesAttractor.} + +\item{theta}{An array of N named dimensions containing the inverse of the +persistence 'theta' as the output of CST_ProxiesAttractor or ProxiesAttractor.} + +\item{ncores}{The number of cores to use in parallel computation} +} +\value{ +A list of length 2: +\itemize{ +\item\code{pred.dim} {a list of two lists 'qdim' and 'pos.d'. The 'qdim' list +contains values of local dimension 'dim' divided by terciles: +d1: lower tercile (high predictability), +d2: middle tercile, +d3: higher tercile (low predictability) +The 'pos.d' list contains the position of each tercile in parameter 'dim'} + +\item\code{pred.theta} {a list of two lists 'qtheta' and 'pos.t'. +The 'qtheta' list contains values of the inverse of persistence 'theta' +divided by terciles: +th1: lower tercile (high predictability), +th2: middle tercile, +th3: higher tercile (low predictability) +The 'pos.t' list contains the position of each tercile in parameter 'theta'} +} + +dyn_scores values from 0 to 1. A dyn_score of 1 indicates the highest +predictability. +} +\description{ +This function divides in terciles the two dynamical proxies +computed with CST_ProxiesAttractor or ProxiesAttractor. These terciles will +be used to measure the predictability of the system in dyn_scores. When the +local dimension 'dim' is small and the inverse of persistence 'theta' is +small the predictability is high, and viceversa. +} +\examples{ +# Creating an example of matrix dat(time,grids): +m <- matrix(rnorm(2000) * 10, nrow = 50, ncol = 40) +names(dim(m)) <- c('time', 'grid') +# imposing a threshold + quanti <- 0.90 +# computing dyn_scores from parameters dim and theta of the attractor +attractor <- ProxiesAttractor(dat = m, quanti = 0.60) +predyn <- Predictability(dim = attractor$dim, theta = attractor$theta) +} +\references{ +Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., +and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large +scale atmospheric predictability.Nature Communications, 10(1), 1316. +DOI = https://doi.org/10.1038/s41467-019-09305-8 " + +Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +Dynamical proxies of North Atlantic predictability and extremes. +Scientific Reports, 7-41278, 2017. +} +\author{ +Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} + +Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} + +Veronica Torralba, \email{veronica.torralba@cmcc.it} + +Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +} diff --git a/man/ProxiesAttractor.Rd b/man/ProxiesAttractor.Rd new file mode 100644 index 0000000000000000000000000000000000000000..768ba7366482ed16390158a991bce22998cf0385 --- /dev/null +++ b/man/ProxiesAttractor.Rd @@ -0,0 +1,61 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_ProxiesAttractor.R +\name{ProxiesAttractor} +\alias{ProxiesAttractor} +\title{Computing two dinamical proxies of the attractor.} +\usage{ +ProxiesAttractor(data, quanti, ncores = NULL) +} +\arguments{ +\item{data}{a multidimensional array with named dimensions to create the attractor. It requires a temporal dimension named 'time' and spatial dimensions called 'lat' and 'lon', or 'latitude' and 'longitude' or 'grid'.} + +\item{quanti}{a number lower than 1 indicating the quantile to perform the computation of local dimension and theta} + +\item{ncores}{The number of cores to use in parallel computation.} +} +\value{ +dim and theta +} +\description{ +This function computes two dinamical proxies of the attractor: +The local dimension (d) and the inverse of the persistence (theta). +These two parameters will be used as a condition for the computation of dynamical +scores to measure predictability and to compute bias correction conditioned by +the dynamics with the function DynBiasCorrection. +Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in: +} +\examples{ +# Example 1: Computing the attractor using simple data +# Creating an example of matrix data(time,grids): +mat <- array(rnorm(36 * 40), c(time = 36, grid = 40)) +qm <- 0.90 # imposing a threshold +Attractor <- ProxiesAttractor(data = mat, quanti = qm) +# to plot the result +time = c(1:length(Attractor$theta)) +layout(matrix(c(1, 3, 2, 3), 2, 2)) +plot(time, Attractor$dim, xlab = 'time', ylab = 'd', + main = 'local dimension', type = 'l') +plot(time, Attractor$theta, xlab = 'time', ylab = 'theta', main = 'theta') +plot(Attractor$dim, Attractor$theta, col = 'blue', + main = "Proxies of the Attractor", + xlab = "local dimension", ylab = "theta", lwd = 8, 'p') + +} +\references{ +Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019). +The hammam effect or how a warm ocean enhances large scale atmospheric predictability. +Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 " + +Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +Dynamical proxies of North Atlantic predictability and extremes. +Scientific Reports, 7-41278, 2017. +} +\author{ +Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} + +Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} + +Veronica Torralba, \email{veronica.torralba@cmcc.it} + +Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +} diff --git a/man/SaveExp.Rd b/man/SaveExp.Rd index 40ace2dbab6cc9bb1c84292ce7387f6a3e991dc0..345974b22a80c52f5cce724fd4f5c560ea876834 100644 --- a/man/SaveExp.Rd +++ b/man/SaveExp.Rd @@ -15,7 +15,8 @@ SaveExp( Dates, cdo_grid_name, projection, - destination + destination, + extra_string = NULL ) } \arguments{ @@ -40,9 +41,11 @@ SaveExp( \item{projection}{a character string indicating the projection name} \item{destination}{a character string indicating the path where to store the NetCDF files} + +\item{extra_string}{a character string to be include as part of the file name, for instance, to identify member or realization.} } \value{ -the function creates as many files as sdates per dataset. Each file could contain multiple members +the function creates as many files as sdates per dataset. Each file could contain multiple members. It would be added to the file name between underscore characters. The path will be created with the name of the variable and each Datasets. } \description{ diff --git a/tests/testthat/test-CST_BiasCorrection.R b/tests/testthat/test-CST_BiasCorrection.R index c5525c57012ec0b92c46e7bc9e2284831a26a66e..fc6405ce7f9ee002c70aa7d27fbb8417bd854e35 100644 --- a/tests/testthat/test-CST_BiasCorrection.R +++ b/tests/testthat/test-CST_BiasCorrection.R @@ -42,4 +42,9 @@ test_that("Sanity checks", { expect_warning(CST_BiasCorrection(exp = exp2, obs = obs2), "Parameter 'obs' contains NA values", "Parameter 'exp' contains NA values.") + + hinc <- array(1:6, c(sdate = 3, member = 2)) + obs <- array(3:6, c(sdate = 3, member = 1)) + res <- round(BiasCorrection(exp = hinc, obs = obs, exp_cor = hinc), 2) + expect_equal(res, array(c(2.66, 4.27, 3.2, 4.8, 3.73, 5.34), c(member = 2, sdate = 3))) }) diff --git a/vignettes/Analogs_vignette.Rmd b/vignettes/Analogs_vignette.Rmd index 5a52a05d99666c617ccf02a1a953b8c60e79c540..31fb9c29842942672502356b7083c6b910922480 100644 --- a/vignettes/Analogs_vignette.Rmd +++ b/vignettes/Analogs_vignette.Rmd @@ -156,7 +156,7 @@ dateseq <- format(seq(start, end, by = "year"), "%Y%m%d") Using the `CST_Load` function from **CSTool package**, the data available in our data store can be loaded. The following lines show how this function can be used. The experimental datasets are interpolated to the ERA5 grid by specifying the 'grid' parameter while ERA5 doesn't need to be interpolated. While parameter leadtimemax is set to 1 for the experimental dataset, it is set to 31 for the observations, returning the daily observations for October for the years requested in 'sdate' (2000-2006). -Download the data to run the recipe in the link https://downloads.cmcc.bo.it/d_chaves/ANALOGS/data_for_Analogs.Rdat or ask carmen.alvarez-castro at cmcc.it or nuria.perez at bsc.es. +Download the data to run the recipe under the HTTPS: downloads.cmcc.bo.it/d_chaves/ANALOGS/data_for_Analogs.Rdat or ask carmen.alvarez-castro at cmcc.it or nuria.perez at bsc.es. ``` exp <- list(name = 'ECMWF_system4_m1', @@ -373,4 +373,4 @@ down4 <- CST_Analogs(expL = expPSL, obsL = obsPSL, AnalogsInfo = TRUE, In this case, the best analog is still being 7th of October, 2005. -*Note: You can compute the anomalies values before applying the criterias (as in Yiou et al, 2013) using `CST_Anomaly` of CSTools package* \ No newline at end of file +*Note: You can compute the anomalies values before applying the criterias (as in Yiou et al, 2013) using `CST_Anomaly` of CSTools package* diff --git a/vignettes/RainFARM_vignette.Rmd b/vignettes/RainFARM_vignette.Rmd index dbcb48a47bee31c92b80c4b879d55d989132843b..c47d0e73b8f03936345e0da19c38630fc0badf06 100644 --- a/vignettes/RainFARM_vignette.Rmd +++ b/vignettes/RainFARM_vignette.Rmd @@ -118,7 +118,7 @@ RainFARM has downscaled the original field with a realistic fine-scale correlati The area of interest in our example presents a complex orography, but the basic RainFARM algorithm used does not consider topographic elevation in deciding how to distribute fine-scale precipitation. A long term climatology of the downscaled fields would have a resolution comparable to that of the original coarse fields and would not resemble the fine-scale structure of an observed climatology. If an external fine-scale climatology of precipitation is available, we can use the method discussed in Terzago et al. (2018) to change the distribution of precipitation by RainFARM for each timestep, so that the long-term average is close to this reference climatology in terms of precipitation distribution (while the total precipitation amount of the original fields to downscale is preserved). -Suitable climatology files could be for example a fine-scale precipitation climatology from a high-resolution regional climate model (see e.g. Terzago et al. 2018), a local high-resolution gridded climatology from observations, or a reconstruction such as those which can be downloaded from the WORLDCLIM (https://www.worldclim.org) or CHELSA (https://chelsa-climate.org) websites. The latter data will need to be converted to NetCDF format before being used (see for example the GDAL tools (https://gdal.org). +Suitable climatology files could be for example a fine-scale precipitation climatology from a high-resolution regional climate model (see e.g. Terzago et al. 2018), a local high-resolution gridded climatology from observations, or a reconstruction such as those which can be downloaded from the WORLDCLIM (https://www.worldclim.org) or CHELSA (chelsa-climate.org) websites. The latter data will need to be converted to NetCDF format before being used (see for example the GDAL tools (https://gdal.org). We will assume that a copy of the WORLDCLIM precipitation climatology at 30 arcseconds (about 1km resolution) is available in the local file `medscope.nc`. From this file we can derive suitable weights to be used with RainFARM using the `CST_RFWeights` functions as follows: ```{r} ww <- CST_RFWeights("./worldclim.nc", nf = 20, lon = exp$lon, lat = exp$lat)