Title:	Assessing Skill of Climate Forecasts on Seasonal-to-Decadal Timescales
Version:	5.3.0
Description:	Exploits dynamical seasonal forecasts in order to provide information relevant to stakeholders at the seasonal timescale. The package contains process-based methods for forecast calibration, bias correction, statistical and stochastic downscaling, optimal forecast combination and multivariate verification, as well as basic and advanced tools to obtain tailored products. This package was developed in the context of the ERA4CS project MEDSCOPE and the H2020 S2S4E project and includes contributions from ArticXchange project founded by EU-PolarNet 2. Implements methods described in Pérez-Zanón et al. (2022) <doi:10.5194/gmd-15-6115-2022>, Doblas-Reyes et al. (2005) <doi:10.1111/j.1600-0870.2005.00104.x>, Mishra et al. (2018) <doi:10.1007/s00382-018-4404-z>, Sanchez-Garcia et al. (2019) <doi:10.5194/asr-16-165-2019>, Straus et al. (2007) <doi:10.1175/JCLI4070.1>, Terzago et al. (2018) <doi:10.5194/nhess-18-2825-2018>, Torralba et al. (2017) <doi:10.1175/JAMC-D-16-0204.1>, D'Onofrio et al. (2014) <doi:10.1175/JHM-D-13-096.1>, Verfaillie et al. (2017) <doi:10.5194/gmd-10-4257-2017>, Van Schaeybroeck et al. (2019) <doi:10.1016/B978-0-12-812372-0.00010-8>, Yiou et al. (2013) <doi:10.1007/s00382-012-1626-3>.
Depends:	R (≥ 3.5.0), maps, qmap, easyVerification
Imports:	s2dv, startR, rainfarmr, multiApply (≥ 2.1.1), ClimProjDiags, ncdf4, plyr, abind, data.table, reshape2, ggplot2, RColorBrewer, graphics, grDevices, stats, utils, verification, lubridate, scales, easyNCDF, dplyr
Suggests:	zeallot, testthat, knitr, markdown, rmarkdown
VignetteBuilder:	knitr
License:	GPL-3
Encoding:	UTF-8
LazyData:	true
RoxygenNote:	7.3.1
Config/testthat/edition:	3
NeedsCompilation:	yes
Packaged:	2025-11-14 10:24:46 UTC; tkariyat
Author:	Nuria Perez-Zanon [aut], Louis-Philippe Caron [aut], Carmen Alvarez-Castro [aut], Lauriane Batte [aut], Carlos Delgado [aut], Jost von Hardenberg [aut], Llorenç LLedo [aut], Nicolau Manubens [aut], Lluís Palma [aut], Eroteida Sanchez-Garcia [aut], Bert van Schaeybroeck [aut], Veronica Torralba [aut], Deborah Verfaillie [aut], Eva Rifà [ctb], Filippo Cali Quaglia [ctb], Maria M. Chaves-Montero [ctb], Chihchung Chou [ctb], Nicola Cortesi [ctb], Susanna Corti [ctb], Paolo Davini [ctb], Gildas Dayon [ctb], Marta Dominguez [ctb], Federico Fabiano [ctb], Ignazio Giuntoli [ctb], Raul Marcos [ctb], Paola Marson [ctb], Niti Mishra [ctb], Jesus Peña [ctb], Francesc Roura-Adserias [ctb], Silvia Terzago [ctb], Danila Volpi [ctb], An-Chi Ho [ctb], Victoria Agudetse [ctb], Theertha Kariyathan [cre], BSC-CNS [cph]
Maintainer:	Theertha Kariyathan <theertha.kariyathan@bsc.es>
Repository:	CRAN
Date/Publication:	2025-11-14 12:40:27 UTC

AdamontQQCorr computes quantile-quantile correction of seasonal or decadal forecast data using weather types

Description

This function computes a quantile mapping based on weather types for experiment data (typically a hindcast) onto reference obs, typically provided by reanalysis data.

Usage

AdamontQQCorr(
  exp,
  wt_exp,
  obs,
  wt_obs,
  corrdims = c("member", "sdate", "ftime"),
  londim = "lon",
  latdim = "lat",
  regrid = FALSE,
  NN = NULL
)

Arguments

Value

An array (such as $data array from an object of class s2dv_cube) with named dimensions, containing experiment data on the lat/lon grid of obs array, corrected by quantile mapping depending on the weather types wt_exp

Author(s)

Paola Marson, paola.marson@meteo.fr for PROSNOW version

Lauriane Batté, lauriane.batte@meteo.fr for CSTools adaptation

Examples

wt_exp <- c(1,1,2,3,3,2,2,1,1,2,2,3)
dim(wt_exp) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3)
wt_obs <- c(3,3,1,2,2,2,2,1,3,1,1,2) 
dim(wt_obs) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3)
exp <- 1 : c(1 * 1 * 4 * 3 * 4 * 4)
dim(exp) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3,
             lat = 4, lon = 4)
obs <- 101 : c(100 + 1 * 1 * 4 * 3 * 4 * 4)
dim(obs) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3,
             lat = 4, lon = 4)
exp_corr <- AdamontQQCorr(exp = exp, wt_exp = wt_exp, 
                         obs = obs, wt_obs = wt_obs, 
                         corrdims = c('dataset', 'member', 'sdate', 'ftime'))

Analogs based on large scale fields.

Description

This function perform a downscaling using Analogs. To compute the analogs, the function search for days with similar large scale conditions to downscaled fields in the local scale. The large scale and the local scale regions are defined by the user. The large scale is usually given by atmospheric circulation as sea level pressure or geopotential height (Yiou et al, 2013) but the function gives the possibility to use another field. The local scale will be usually given by precipitation or temperature fields, but might be another variable. The analogs function will find the best analogs based in three criterias: (1) Minimum Euclidean distance in the large scale pattern (i.e. SLP) (2) Minimum Euclidean distance in the large scale pattern (i.e. SLP) and minimum Euclidean distance in the local scale pattern (i.e. SLP). (3) Minimum Euclidean distance in the large scale pattern (i.e. SLP), minimum distance in the local scale pattern (i.e. SLP) and highest correlation in the local variable to downscale (i.e Precipitation). The search of analogs must be done in the longest dataset posible. This is important since it is necessary to have a good representation of the possible states of the field in the past, and therefore, to get better analogs. Once the search of the analogs is complete, and in order to used the three criterias the user can select a number of analogs , using parameter 'nAnalogs' to restrict the selection of the best analogs in a short number of posibilities, the best ones. This function has not constrains of specific regions, variables to downscale, or data to be used (seasonal forecast data, climate projections data, reanalyses data). The regrid into a finner scale is done interpolating with CST_Start. Then, this interpolation is corrected selecting the analogs in the large and local scale in based of the observations. The function is an adapted version of the method of Yiou et al 2013.

Usage

Analogs(
  expL,
  obsL,
  time_obsL,
  time_expL = NULL,
  lonL = NULL,
  latL = NULL,
  expVar = NULL,
  obsVar = NULL,
  sdate_dim = "sdate",
  criteria = "Large_dist",
  excludeTime = NULL,
  lonVar = NULL,
  latVar = NULL,
  region = NULL,
  nAnalogs = NULL,
  AnalogsInfo = FALSE,
  ncores = NULL
)

Arguments

Value

An array with the dowscaled values of the best analogs for the criteria selected. If 'AnalogsInfo' is set to TRUE it returns a list with an array of the dowsncaled field and the analogs information in elements 'analogs', 'metric' and 'dates'.

Author(s)

M. Carmen Alvarez-Castro, carmen.alvarez-castro@cmcc.it

Maria M. Chaves-Montero, mariadm.chaves@cmcc.it

Veronica Torralba, veronica.torralba@cmcc.it

Nuria Perez-Zanon nuria.perez@bsc.es

References

Yiou, P., T. Salameh, P. Drobinski, L. Menut, R. Vautard, and M. Vrac, 2013 : Ensemble reconstruction of the atmospheric column from surface pressure using analogues. Clim. Dyn., 41, 1419-1437. pascal.yiou@lsce.ipsl.fr

Examples

# Example 1: Downscaling using criteria 'Large_dist' and a single variable:
expSLP <- rnorm(1:20)
dim(expSLP) <- c(lat = 4, lon = 5)
obsSLP <- c(rnorm(1:180), expSLP * 1.2)
dim(obsSLP) <- c(time = 10, lat = 4, lon = 5)
time_obsSLP <- paste(rep("01", 10), rep("01", 10), 1994 : 2003, sep = "-")
dim(time_obsSLP) <- c(time = 10)
downscale_field <- Analogs(expL = expSLP, obsL = obsSLP, 
                          time_obsL = time_obsSLP,time_expL = "01-01-1994")

# Example 2: Downscaling using criteria 'Large_dist' and 2 variables:
obs.pr <- c(rnorm(1:200) * 0.001)
dim(obs.pr) <- dim(obsSLP)
downscale_field <- Analogs(expL = expSLP, obsL = obsSLP, obsVar = obs.pr,
                          time_obsL = time_obsSLP, time_expL = "01-01-1994")

# Example 3: Downscaling using criteria 'Local_dist' and 2 variables:
# analogs of local scale using criteria 2
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", 
                      lonL = seq(-1, 5, 1.5),latL = seq(30, 35, 1.5), 
                      region = region,time_expL = "01-10-2000", 
                      nAnalogs = 10, AnalogsInfo = TRUE)

# Example 4: Downscaling using criteria 'Local_cor' and 2 variables:
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", 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)

# Example 5: List of best analogs in the three criterias Large_dist, 
Large_scale <- Analogs(expL = expSLP, obsL = obsSLP, time_obsL = time_obsSLP,
                      criteria = "Large_dist", time_expL = "01-10-2000",
                      nAnalogs = 7, AnalogsInfo = TRUE)
Local_scale <- Analogs(expL = expSLP, obsL = obsSLP, time_obsL = time_obsSLP,
                      time_expL = "01-10-2000", criteria = "Local_dist",
                      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)

Calculate the spatial area-weighted average of multidimensional arrays.

Description

This function computes a spatial area-weighted average of n-dimensional arrays given a spatial mask containing the area (e.g. in m^2) for each grid box.

Usage

AreaWeighted(data, area, lon_dim = "lon", lat_dim = "lat", ncores = NULL)

Arguments

Value

An array with the same dimensions as data, except with lon_dim and lat_dim removed. Any extra dimensions in area not present in data are preserved in the output (e.g., region dimension).

Examples

data <- array(data = 1:10, dim = c(year = 10, lon = 5, lat = 3))
area <- array(data = 1:10, dim = c(lon = 5, lat = 3))
AreaWeighted(data, area)
# With extra region dimension
area <- array(data = 1:10, dim = c(lon = 5, lat = 3, region = 4))
AreaWeighted(data, area)

Computing the weighted ensemble means for SFSs.

Description

This function implements the computation to obtain the weighted ensemble means for SFSs using a normalized weights array,

Usage

BEI_EMWeighting(var_exp, aweights, time_dim_name = "time", memb_dim = "member")

Arguments

Value

BEI_EMWeighting() returns an array with at least one or three dimensions depending if the variable is spatially aggregated variable (as e.g. NAO index)(time) or it is spatial variable (as e.g. precipitation or temperature) (time, lat, lon), containing the ensemble means computing with weighted members.

Author(s)

Eroteida Sanchez-Garcia - AEMET, esanchezg@aemet.es

References

Regionally improved seasonal forecast of precipitation through Best estimation of winter NAO, Sanchez-Garcia, E. et al., Adv. Sci. Res., 16, 165174, 2019, doi:10.5194/asr-16-165-2019

Examples

# Example 1 
var_exp <- 1 : (2 * 3 * 4)
dim(var_exp) <- c(time = 2, dataset = 3, member = 4)
aweights <- runif(24, min = 0.001, max = 0.999)
dim(aweights) <- c(time = 2, dataset = 3, member = 4)
res <- BEI_EMWeighting(var_exp, aweights)

# Example 2 
var_exp <- 1 : (2 * 4 * 2 * 3)
dim(var_exp) <- c(time = 2, member = 4, lat = 2, lon = 3)
aweights <- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3)
dim(aweights) <- c(time = 2, member = 4)
res <- BEI_EMWeighting(var_exp, aweights)

Computing the Best Index PDFs combining Index PDFs from two SFSs

Description

This function implements the computation to obtain the index Probability Density Functions (PDFs) (e.g. NAO index) obtained to combining the Index PDFs for two Seasonal Forecast Systems (SFSs), the Best Index estimation (see Sanchez-Garcia, E. et al (2019), doi:10.5194/asr-16-165-2019 for more details about the methodology applied to estimate the Best Index).

Usage

BEI_PDFBest(
  index_obs,
  index_hind1,
  index_hind2,
  index_fcst1 = NULL,
  index_fcst2 = NULL,
  method_BC = "none",
  time_dim_name = "time",
  na.rm = FALSE
)

Arguments

Value

BEI_PDFBest() returns an array with the parameters that caracterize the PDFs, with at least a temporal dimension, by default 'time' and dimension 'statistic' equal to 2. The firt statistic is the parameter 'mean' of the PDF for the best estimation combining the two SFSs PDFs. The second statistic is the parameter 'standard deviation' of the PDF for the best estimation combining the two SFSs PDFs. If index_fcst1 and/or index_fcst2 are null, returns the values for hindcast period. Otherwise, it returns the values for a forecast year.

Author(s)

Eroteida Sanchez-Garcia - AEMET, esanchezg@aemet.es

References

Regionally improved seasonal forecast of precipitation through Best estimation of winter NAO, Sanchez-Garcia, E. et al., Adv. Sci. Res., 16, 165174, 2019, doi:10.5194/asr-16-165-2019

Examples

# Example 1 for the BEI_PDFBest function
index_obs<- rnorm(10, sd = 3)
dim(index_obs) <- c(time = 5, season = 2)
index_hind1 <- rnorm(40, mean = 0.2, sd = 3)
dim(index_hind1) <- c(time = 5, member = 4, season = 2)
index_hind2 <- rnorm(60, mean = -0.5, sd = 4)
dim(index_hind2) <- c(time = 5, member = 6, season = 2)
index_fcst1 <- rnorm(16, mean = 0.2, sd = 3)
dim(index_fcst1) <- c(time = 1, member = 8, season = 2)
index_fcst2 <- rnorm(18, mean = -0.5, sd = 4)
dim(index_fcst2) <- c(time = 1, member = 9, season = 2)
method_BC <- 'ME'
res <- BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, 
index_fcst2, method_BC)  
# Example 2 for the BEI_PDFBest function
index_obs<- rnorm(10, sd = 3)
dim(index_obs) <- c(time = 5, season = 2)
index_hind1 <- rnorm(40, mean = 0.2, sd = 3)
dim(index_hind1) <- c(time = 5, member = 4, season = 2)
index_hind2 <- rnorm(60, mean = -0.5, sd = 4)
dim(index_hind2) <- c(time = 5, member = 6, season = 2)
index_fcst1 <- rnorm(16, mean = 0.2, sd = 3)
dim(index_fcst1) <- c(time = 1, member = 8, season = 2)
index_fcst2 <- rnorm(18, mean = -0.5, sd = 4)
dim(index_fcst2) <- c(time = 1, member = 9, season = 2)
method_BC <- c('LMEV', 'ME')
res <- BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, 
                   index_fcst2, method_BC) 

Computing the weighted tercile probabilities for SFSs.

Description

This function implements the computation to obtain the tercile probabilities for a weighted variable for SFSs using a normalized weights array,

Usage

BEI_ProbsWeighting(
  var_exp,
  aweights,
  terciles,
  time_dim_name = "time",
  memb_dim = "member"
)

Arguments

Value

BEI_ProbsWeighting() returns an array with at least two or four dimensions depending if the variable is a spatially aggregated variable (as e.g. NAO index)(time, tercil) or it is spatial variable (as e.g. precipitation or temperature)(time, tercile, lat, lon), containing the terciles probabilities computing with weighted members. The first tercil is the lower tercile, the second is the normal tercile and the third is the upper tercile.

Author(s)

Eroteida Sanchez-Garcia - AEMET, esanchezg@aemet.es

References

Regionally improved seasonal forecast of precipitation through Best estimation of winter NAO, Sanchez-Garcia, E. et al., Adv. Sci. Res., 16, 165174, 2019, doi:10.5194/asr-16-165-2019

Examples

# Example 1 
var_exp <- 1 : (2 * 4)
dim(var_exp) <- c(time = 2, member = 4)
aweights <- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3)
dim(aweights) <- c(time = 2, member = 4)
terciles <- c(2.5,5)
dim(terciles) <- c(tercil = 2)
res <- BEI_ProbsWeighting(var_exp, aweights, terciles)

# Example 2 
var_exp <- rnorm(48, 50, 9)
dim(var_exp) <- c(time = 2, member = 4, lat = 2, lon = 3)
aweights <- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3)
dim(aweights) <- c(time = 2, member = 4)
terciles <- rep(c(48,50), 2*3)
dim(terciles) <- c(tercil = 2, lat = 2, lon = 3)
res <- BEI_ProbsWeighting(var_exp, aweights, terciles)

Computing the weighted terciles for SFSs.

Description

This function implements the computation to obtain the terciles for a weighted variable for SFSs using a normalized weights array,

Usage

BEI_TercilesWeighting(
  var_exp,
  aweights,
  time_dim_name = "time",
  memb_dim = "member"
)

Arguments

Value

BEI_TercilesWeighting() returns an array with at least one dimension depending if the variable is a spatially aggregated variable (as e.g. NAO index)(tercil) or it is spatial variable (as e.g. precipitation or temperature)(tercil, lat, lon), containing the terciles computing with weighted members. The first tercil is the lower tercile, the second is the upper tercile.

Author(s)

Eroteida Sanchez-Garcia - AEMET, esanchezg@aemet.es

References

Regionally improved seasonal forecast of precipitation through Best estimation of winter NAO, Sanchez-Garcia, E. et al., Adv. Sci. Res., 16, 165174, 2019, doi:10.5194/asr-16-165-2019

Examples

# Example 1 
var_exp <- 1 : (2 * 4)
dim(var_exp) <- c(time = 2, member = 4)
aweights<- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3)
dim(aweights) <- c(time = 2, member = 4)
res <- BEI_TercilesWeighting(var_exp, aweights)

# Example 2 
var_exp <- rnorm(48, 50, 9)
dim(var_exp) <- c(time = 2, member = 4, lat = 2, lon = 3)
aweights<- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3)
dim(aweights) <- c(time = 2, member = 4)
res <- BEI_TercilesWeighting(var_exp, aweights)

Computing the weights for SFSs using the Best Index PDFs.

Description

This function implements the computation to obtain the normalized weights for each member of each Seasonal Forecast Systems (SFS) or dataset using the Probability Density Functions (PDFs) indicated by the parameter 'pdf_weight' (for instance the Best Index estimation obtained using the 'PDFBest' function). The weight of each member is proportional to the probability of its index calculated with the PDF "pdf_weight".

Usage

BEI_Weights(index_weight, pdf_weight, time_dim_name = "time")

Arguments

Value

BEI_Weights() returns a normalized weights array with the same dimensions that index_weight.

Author(s)

Eroteida Sanchez-Garcia - AEMET, esanchezg@aemet.es

References

Regionally improved seasonal forecast of precipitation through Best estimation of winter NAO, Sanchez-Garcia, E. et al., Adv. Sci. Res., 16, 165174, 2019, doi:10.5194/asr-16-165-2019

Examples

# Example for the BEI_Weights function
index_weight <- 1 : (10 * 3 * 5 * 1)
dim(index_weight) <- c(sdate = 10, dataset = 3, member = 5, season = 1)
pdf_weight <- 1 : (10 * 3 * 2 * 1)
dim(pdf_weight) <- c(sdate = 10, dataset = 3, statistic = 2, season = 1)
res <- BEI_Weights(index_weight, pdf_weight)
dim(res)
# sdate   dataset    member season
#    10         3         5      1

Bias Correction based on the mean and standard deviation adjustment

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.

Usage

BiasCorrection(
  exp,
  obs,
  exp_cor = NULL,
  na.rm = FALSE,
  memb_dim = "member",
  sdate_dim = "sdate",
  dat_dim = NULL,
  ncores = NULL
)

Arguments

Value

An array containing the bias corrected forecasts with the dimensions nexp, nobs and same dimensions as in the 'exp' object. nexp is the number of experiment (i.e., 'dat_dim' in exp), and nobs is the number of observation (i.e., 'dat_dim' in obs). If dat_dim is NULL, nexp and nobs are omitted. If 'exp_cor' is provided the returned array will be with the same dimensions as 'exp_cor'.

Author(s)

Verónica Torralba, veronica.torralba@bsc.es

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)

Examples

mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7)
dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, time = 5, lat = 6, lon = 7)
obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, time = 5, lat = 6, lon = 7)
a <- BiasCorrection(exp = mod1, obs = obs1)

Bind two arrays by a specified named dimension

Description

This function combines the data inside two or more arrays with named dimensions along the specified along dimension. It is a wrapper of the abind() function from the abind package.

Usage

BindDim(x, along)

Arguments

Value

A single array combining the arrays in x along the specified dimension, with dimension names. The order of the dimensions will be the same as in the first array provided in the list.

Author(s)

Agudetse Roures Victoria, victoria.agudetse@bsc.es

See Also

abind

Examples

array1 <- array(1:100, dim = c(time = 1, lon = 10, lat = 10))
array2 <- array(101:200, dim = c(lon = 10, lat = 10, time = 1))
# Bind arrays
array3 <- BindDim(x = list(array1, array2),
                 along = "time")
# Check new dimensions
dim(array3)

CST_AdamontAnalog finds analogous data in the reference dataset to experiment data based on weather types

Description

This function searches for analogs in a reference dataset for experiment data, based on corresponding weather types. The experiment data is typically a hindcast, observations are typically provided by reanalysis data.

This function searches for analogs in a reference dataset for experiment data, based on corresponding weather types. The experiment data is typically a hindcast, observations are typically provided by reanalysis data.

Usage

CST_AdamontAnalog(
  exp,
  obs,
  wt_exp,
  wt_obs,
  nanalogs,
  method = "pattcorr",
  thres = NULL,
  search_obsdims = c("member", "sdate", "ftime"),
  londim = "lon",
  latdim = "lat"
)

AdamontAnalog(
  exp,
  obs,
  wt_exp,
  wt_obs,
  nanalogs = 5,
  method = "pattcorr",
  thres = NULL,
  search_obsdims = c("member", "sdate", "ftime"),
  londim = "lon",
  latdim = "lat"
)

Arguments

Value

analog_vals An object of class s2dv_cube containing nanalogs analog values for each value of exp input data.

analog_vals An array containing nanalogs analog values.

Author(s)

Paola Marson, paola.marson@meteo.fr for PROSNOW version

Lauriane Batté, lauriane.batte@meteo.fr for CSTools adaptation

Examples

wt_exp <- sample(1:3, 15*6*3, replace = TRUE)
dim(wt_exp) <- c(dataset = 1, member = 15, sdate = 6, ftime = 3)
wt_obs <- sample(1:3, 6*3, replace = TRUE)
dim(wt_obs) <- c(dataset = 1, member = 1, sdate = 6, ftime = 3)
exp <- NULL
exp$data <- 1 : c(1 * 15 * 6 * 3 * 8 * 8)
dim(exp$data) <- c(dataset = 1, member = 15, sdate = 6, ftime = 3,
                  lat = 8, lon = 8)
class(exp) <- 's2dv_cube'
obs <- NULL
obs$data <- 101 : c(100 + 1 * 1 * 6 * 3 * 8 * 8)
dim(obs$data) <- c(dataset = 1, member = 1, sdate = 6, ftime = 3,
                  lat = 8, lon = 8)
class(obs) <- 's2dv_cube'
analog_vals <- CST_AdamontAnalog(exp = exp, obs = obs, wt_exp = wt_exp, 
                                wt_obs = wt_obs, nanalogs = 2)
wt_exp <- sample(1:3, 15*6*3, replace = TRUE)
dim(wt_exp) <- c(dataset = 1, member = 15, sdate = 6, ftime = 3)
wt_obs <- sample(1:3, 6*3, replace = TRUE)
dim(wt_obs) <- c(dataset = 1, member = 1, sdate = 6, ftime = 3)
exp <- 1 : c(1 * 15 * 6 * 3 * 8 * 8)
dim(exp) <- c(dataset = 1, member = 15, sdate = 6, ftime = 3, lat = 8, lon = 8)
obs <- 101 : c(100 + 1 * 1 * 6 * 3 * 8 * 8)
dim(obs) <- c(dataset = 1, member = 1, sdate = 6, ftime = 3, lat = 8, lon = 8)
analog_vals <- AdamontAnalog(exp = exp, obs = obs, wt_exp = wt_exp, 
                            wt_obs = wt_obs, nanalogs = 2)

CST_AdamontQQCorr computes quantile-quantile correction of seasonal or decadal forecast data using weather types

Description

This function computes a quantile mapping based on weather types for experiment data (typically a hindcast) onto reference obs, typically provided by reanalysis data.

Usage

CST_AdamontQQCorr(
  exp,
  wt_exp,
  obs,
  wt_obs,
  corrdims = c("member", "sdate", "ftime"),
  londim = "lon",
  latdim = "lat"
)

Arguments

Value

An object of class s2dv_cube containing experiment data on the lat/lon grid of obs input data, corrected by quantile mapping depending on the weather types wt_exp.

Author(s)

Lauriane Batté, lauriane.batte@meteo.fr

Paola Marson, paola.marson@meteo.fr

Gildas Dayon, gildas.dayon@meteo.fr

Examples

wt_exp <- c(1,1,2,3,3,2,2,1,1,2,2,3)
dim(wt_exp) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3)
wt_obs <- c(3,3,1,2,2,2,2,1,3,1,1,2) 
dim(wt_obs) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3)
exp <- NULL
exp$data <-  1 : c(1 * 1 * 4 * 3 * 4 * 4)
dim(exp$data) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3,
                  lat = 4, lon = 4)
class(exp) <- 's2dv_cube'
obs <- NULL
obs$data <-  101 : c(100 + 1 * 1 * 4 * 3 * 4 * 4)
dim(obs$data) <- c(dataset = 1, member = 1, sdate = 4, ftime = 3,
                  lat = 4, lon = 4)
class(obs) <- 's2dv_cube'
exp_corr <- CST_AdamontQQCorr(exp = exp, wt_exp = wt_exp, 
                             obs = obs, wt_obs = wt_obs, 
                             corrdims = c('dataset','member','sdate','ftime'))

Downscaling using Analogs based on large scale fields.

Description

This function perform a downscaling using Analogs. To compute the analogs, the function search for days with similar large scale conditions to downscaled fields to a local scale. The large scale and the local scale regions are defined by the user. The large scale is usually given by atmospheric circulation as sea level pressure or geopotential height (Yiou et al, 2013) but the function gives the possibility to use another field. The local scale will be usually given by precipitation or temperature fields, but might be another variable.The analogs function will find the best analogs based in Minimum Euclidean distance in the large scale pattern (i.e.SLP).

The search of analogs must be done in the longest dataset posible. This is important since it is necessary to have a good representation of the possible states of the field in the past, and therefore, to get better analogs. This function has not constrains of specific regions, variables to downscale, or data to be used (seasonal forecast data, climate projections data, reanalyses data). The regrid into a finner scale is done interpolating with CST_Start. Then, this interpolation is corrected selecting the analogs in the large and local scale in based of the observations. The function is an adapted version of the method of Yiou et al 2013. For an advanced search of Analogs (multiple Analogs, different criterias, further information from the metrics and date of the selected Analogs) use the'Analog' function within 'CSTools' package.

Usage

CST_Analogs(
  expL,
  obsL,
  expVar = NULL,
  obsVar = NULL,
  sdate_dim = "sdate",
  region = NULL,
  criteria = "Large_dist",
  excludeTime = NULL,
  time_expL = NULL,
  time_obsL = NULL,
  nAnalogs = NULL,
  AnalogsInfo = FALSE,
  ncores = NULL
)

Arguments

Value

An 's2dv_cube' object containing an array with the dowscaled values of the best analogs in element 'data'. If 'AnalogsInfo' is TRUE, 'data' is a list with an array of the downscaled fields and the analogs information in elements 'analogs', 'metric' and 'dates'.

Author(s)

M. Carmen Alvarez-Castro, carmen.alvarez-castro@cmcc.it

Maria M. Chaves-Montero, mariadm.chaves@cmcc.it

Veronica Torralba, veronica.torralba@cmcc.it

Nuria Perez-Zanon nuria.perez@bsc.es

References

Yiou, P., T. Salameh, P. Drobinski, L. Menut, R. Vautard, and M. Vrac, 2013 : Ensemble reconstruction of the atmospheric column from surface pressure using analogues. Clim. Dyn., 41, 1419-1437. pascal.yiou@lsce.ipsl.fr

See Also

CST_Start, Start

Examples

expL <- rnorm(1:200)
dim(expL) <- c(member = 10, lat = 4, lon = 5)
obsL <- c(rnorm(1:180), expL[1, , ]*1.2)
dim(obsL) <- c(time = 10, lat = 4, lon = 5)
time_obsL <- as.POSIXct(paste(rep("01", 10), rep("01", 10), 1994:2003, sep = "-"), 
                       format = "%d-%m-%y")
dim(time_obsL) <- c(time = 10)
time_expL <- time_obsL[1]
dim(time_expL) <- c(time = 1)
lon <-  seq(-1, 5, 1.5)
lat <- seq(30, 35, 1.5)
coords <- list(lon = seq(-1, 5, 1.5), lat = seq(30, 35, 1.5))
attrs_expL <- list(Dates = time_expL)
attrs_obsL <- list(Dates = time_obsL)
expL <- list(data = expL, coords = coords, attrs = attrs_expL)
obsL <- list(data = obsL, coords = coords, attrs = attrs_obsL)
class(expL) <- 's2dv_cube'
class(obsL) <- 's2dv_cube'
region <- c(min(lon), max(lon), min(lat), max(lat))
downscaled_field <- CST_Analogs(expL = expL, obsL = obsL, region = region)

AEMET Downscaling Precipitation and maximum and minimum temperature downscaling method based on analogs: synoptic situations and significant predictors.

Description

This function downscales low resolution precipitation data (e.g. from Seasonal Forecast Models) through the association with an observational high resolution (HR) dataset (AEMET 5 km gridded data of daily precipitation (Peral et al., 2017)) and a collection of predictors and past synoptic situations similar to estimated day. The method uses three domains:

• Peninsular Spain and Balearic Islands domain (5 km resolution): HR precipitation and the downscaling result domain.

• Synoptic domain (low resolution, e.g. 1.5º x 1.5º): it should be centered over Iberian Peninsula and cover enough extension to detect as much synoptic situations as possible.

• Extended domain (low resolution, e.g. 1.5º x 1.5º): it should have the same resolution as synoptic domain. It is used for SLP Seasonal Forecast Models.

Usage

CST_AnalogsPredictors(
  exp,
  slp,
  obs,
  lon,
  lat,
  slp_lon,
  slp_lat,
  var_name,
  hr_obs,
  tdates,
  ddates,
  restrain,
  dim_name_longitude = "lon",
  dim_name_latitude = "lat",
  dim_name_time = "time"
)

Arguments

Value

Matrix with seasonal forecast precipitation (mm) or maximum and minimum temperature (dozens of ºC) in a 5km x 5km regular grid over peninsular Spain and Balearic Islands. The resulted matrices have two dimensions ('ddates' x 'nptos').(ddates = number of downscaling days and nptos = number of 'hr_obs' gridpoints).

Author(s)

Marta Dominguez Alonso - AEMET, mdomingueza@aemet.es

Nuria Perez-Zanon - BSC, nuria.perez@bsc.es

Anomalies relative to a climatology along selected dimension with or without cross-validation

Description

This function computes the anomalies relative to a climatology computed along the selected dimension (usually starting dates or forecast time) allowing the application or not of crossvalidated climatologies. The computation is carried out independently for experimental and observational data products.

Usage

CST_Anomaly(
  exp = NULL,
  obs = NULL,
  dim_anom = "sdate",
  cross = FALSE,
  memb_dim = "member",
  memb = TRUE,
  dat_dim = c("dataset", "member"),
  filter_span = NULL,
  ftime_dim = "ftime",
  ncores = NULL
)

Arguments

Value

A list with two S3 objects, 'exp' and 'obs', of the class 's2dv_cube', containing experimental and date-corresponding observational anomalies, respectively. These 's2dv_cube's can be ingested by other functions in CSTools.

Author(s)

Perez-Zanon Nuria, nuria.perez@bsc.es

Pena Jesus, jesus.pena@bsc.es

See Also

Ano_CrossValid, Clim and CST_Start

Examples

mod <- 1 : (2 * 3 * 4 * 5 * 6 * 7)
dim(mod) <- c(dataset = 2, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
obs <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
dim(obs) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
lon <- seq(0, 30, 5)
lat <- seq(0, 25, 5)
coords <- list(lon = lon, lat = lat)
exp <- list(data = mod, coords = coords)
obs <- list(data = obs, coords = coords)
attr(exp, 'class') <- 's2dv_cube'
attr(obs, 'class') <- 's2dv_cube'

anom <- CST_Anomaly(exp = exp, obs = obs, cross = FALSE, memb = TRUE)

Calculate the spatial area-weighted average of multidimensional arrays.

Description

This function computes a spatial area-weighted average of n-dimensional arrays given a spatial mask containing the area (e.g. in m^2) for each grid box.

Usage

CST_AreaWeighted(
  data,
  area,
  lon_dim = "lon",
  lat_dim = "lat",
  ncores = NULL,
  extra_info = NULL
)

Arguments

Value

An object of class s2dv_cube with same dimensions as in object data, except with lon_dim and lat_dim removed. Any extra dimensions in area not present in data are preserved in the output (e.g., region dimension).

Examples

data <- array(data = 1:10, dim = c(year = 10, lon = 5, lat = 3))
coords <- list(lat = 1:3, lon = 1:5)
data <- list(data = data, coords = coords)
attr(data, 'class') <- 's2dv_cube'
area <- array(data = 1:10, dim = c(lon = 5, lat = 3))
CST_AreaWeighted(data, area)
# With extra region dimension
area <- array(data = 1:10, dim = c(lon = 5, lat = 3, region = 4))
CST_AreaWeighted(data, area)

Weighting SFSs of a CSTools object.

Description

Function to apply weights to a 's2dv_cube' object. It could return a weighted ensemble mean (deterministic output) or the terciles probabilities (probabilistic output) for Seasonal Forecast Systems (SFSs).

Usage

CST_BEI_Weighting(
  var_exp,
  aweights,
  terciles = NULL,
  type = "ensembleMean",
  time_dim_name = "time",
  memb_dim = "member"
)

Arguments

Value

CST_BEI_Weighting() returns a CSTools object (i.e., of the class 's2dv_cube'). This object has at least an element named $data with at least a temporal dimension (and dimension 'tercil' when the output are tercile probabilities), containing the ensemble means computing with weighted members or probabilities of terciles.

Author(s)

Eroteida Sanchez-Garcia - AEMET, esanchezg@aemet.es

References

Regionally improved seasonal forecast of precipitation through Best estimation of winter NAO, Sanchez-Garcia, E. et al., Adv. Sci. Res., 16, 165174, 2019, doi:10.5194/asr-16-165-2019

Examples

var_exp <- 1 : (2 * 4 * 3 * 2)
dim(var_exp) <- c(time = 2, member = 4, lat = 3, lon = 2)
aweights <- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3, 0.1, 0.2, 0.4, 0.4, 0.1, 
             0.2, 0.4, 0.2)
dim(aweights) <- c(time = 2, member = 4, dataset = 2)
var_exp <- list(data = var_exp)
class(var_exp) <- 's2dv_cube'
res_CST <- CST_BEI_Weighting(var_exp, aweights)

Bias Correction based on the mean and standard deviation adjustment

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.

Usage

CST_BiasCorrection(
  exp,
  obs,
  exp_cor = NULL,
  na.rm = FALSE,
  memb_dim = "member",
  sdate_dim = "sdate",
  dat_dim = NULL,
  ncores = NULL
)

Arguments

Value

An object of class s2dv_cube containing the bias corrected forecasts with the dimensions nexp, nobs and same dimensions as in the 'exp' object. nexp is the number of experiment (i.e., 'dat_dim' in exp), and nobs is the number of observation (i.e., 'dat_dim' in obs). If dat_dim is NULL, nexp and nobs are omitted. If 'exp_cor' is provided the returned array will be with the same dimensions as 'exp_cor'.

Author(s)

Verónica Torralba, veronica.torralba@bsc.es

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)

Examples

mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7)
dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, time = 5, lat = 6, lon = 7)
obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, time = 5, lat = 6, lon = 7)
lon <- seq(0, 30, 5)
lat <- seq(0, 25, 5)
coords <- list(lat = lat, lon = lon)
exp <- list(data = mod1, coords = coords)
obs <- list(data = obs1, coords = coords)
attr(exp, 'class') <- 's2dv_cube'
attr(obs, 'class') <- 's2dv_cube'
a <- CST_BiasCorrection(exp = exp, obs = obs)

Bind two objects of class s2dv_cube

Description

This function combines the data inside two or more objects of class s2dv_cube along the specified along dimension, and modifies the dimensions, coordinates and attributes accordingly, producing a result that contains the complete metadata for all variables, time steps and spatial coordinates that are bound in the process. It ensures that the information inside the s2dv_cube remains coherent with the data it contains.

If the dimension specified in along is among the time dimensions in attrs$Dates, the dates arrays are also bound along this dimension. The load_parameters and when attributes of the first cube are preserved. The source_files attribute is bound along the var_dim and dat_dim dimensions.

Usage

CST_BindDim(x, along, var_dim = NULL, dat_dim = NULL)

Arguments

Value

An object of class s2dv_cube with the combined data, dimensions, coordinates and attributes of the elements in x.

Author(s)

Agudetse Roures Victoria, victoria.agudetse@bsc.es

See Also

abind

Examples

 # Example with sample data:
 # Check original dimensions and coordinates
 lonlat_temp$exp$dims
 lonlat_temp$obs$dims
 # Bind both datasets along the member dimension
 res <- CST_BindDim(x = list(lonlat_temp$exp, lonlat_temp$obs),
                    along = "member",
                    var_dim = NULL,
                    dat_dim = "dat")
 # Check new dimensions and coordinates
 res$dims
 names(res$coords)

Forecast Calibration

Description

Five types of member-by-member bias correction can be performed. The "bias" method corrects the bias only, the "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 "mse_min" and "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 "mse_min" method minimizes a constrained mean-squared error using three parameters, the "crps_min" method features four parameters and minimizes the Continuous Ranked Probability Score (CRPS). The "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). It is equivalent to function Calibration but for objects of class s2dv_cube.

Usage

CST_Calibration(
  exp,
  obs,
  exp_cor = NULL,
  cal.method = "mse_min",
  eval.method = "leave-k-out",
  multi.model = FALSE,
  na.fill = TRUE,
  na.rm = TRUE,
  apply_to = NULL,
  alpha = NULL,
  memb_dim = "member",
  sdate_dim = "sdate",
  dat_dim = NULL,
  ncores = NULL,
  k = 1,
  tail.out = TRUE
)

Arguments

Details

Both the na.fill and na.rm parameters can be used to indicate how the function has to handle the NA values. The na.fill parameter checks whether there are more than three forecast-observations pairs to perform the computation. In case there are three or less pairs, the computation is not carried out, and the value returned by the function depends on the value of this parameter (either NA if na.fill == TRUE or the uncorrected value if na.fill == TRUE). On the other hand, na.rm is used to indicate the function whether to remove the missing values during the computation of the parameters needed to perform the calibration.

Value

An object of class s2dv_cube containing the calibrated forecasts in the element data with the dimensions nexp, nobs and same dimensions as in the 'exp' object. nexp is the number of experiment (i.e., 'dat_dim' in exp), and nobs is the number of observation (i.e., 'dat_dim' in obs). If dat_dim is NULL, nexp and nobs are omitted. If 'exp_cor' is provided the returned array will be with the same dimensions as 'exp_cor'.

Author(s)

Verónica Torralba, veronica.torralba@bsc.es

Bert Van Schaeybroeck, bertvs@meteo.be

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

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

Van Schaeybroeck, B., & Vannitsem, S. (2011). Post-processing through linear regression. Nonlinear Processes in Geophysics, 18(2), 147. doi:10.5194/npg-18-147-2011

Van Schaeybroeck, B., & Vannitsem, S. (2015). Ensemble post-processing using member-by-member approaches: theoretical aspects. Quarterly Journal of the Royal Meteorological Society, 141(688), 807-818. doi:10.1002/qj.2397

See Also

CST_Start

Examples

# Example 1:
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)
lon <- seq(0, 30, 5)
lat <- seq(0, 25, 5)
coords <- list(lat = lat, lon = lon)
exp <- list(data = mod1, coords = coords)
obs <- list(data = obs1, coords = coords)
attr(exp, 'class') <- 's2dv_cube'
attr(obs, 'class') <- 's2dv_cube'
a <- CST_Calibration(exp = exp, obs = obs, cal.method = "mse_min", eval.method = "in-sample")

# Example 2:
mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7)
mod2 <- 1 : (1 * 3 * 1 * 5 * 6 * 7)
dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
dim(mod2) <- c(dataset = 1, member = 3, sdate = 1, 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)
lon <- seq(0, 30, 5)
lat <- seq(0, 25, 5)
coords <- list(lat = lat, lon = lon)
exp <- list(data = mod1, coords = coords)
obs <- list(data = obs1, coords = coords)
exp_cor <- list(data = mod2, lat = lat, lon = lon)
attr(exp, 'class') <- 's2dv_cube'
attr(obs, 'class') <- 's2dv_cube'
attr(exp_cor, 'class') <- 's2dv_cube'
a <- CST_Calibration(exp = exp, obs = obs, exp_cor = exp_cor, cal.method = "evmos")

Make categorical forecast based on a multi-model forecast with potential for calibrate

Description

This function converts a multi-model ensemble forecast into a categorical forecast by giving the probability for each category. Different methods are available to combine the different ensemble forecasting models into probabilistic categorical forecasts.

Motivation: Beyond the short range, the unpredictable component of weather predictions becomes substantial due to the chaotic nature of the earth system. Therefore, predictions can mostly be skillful when used in a probabilistic sense. In practice this is done using ensemble forecasts. It is then common to convert the ensemble forecasts to occurence probabilities for different categories. These categories typically are taken as terciles from climatolgical distributions. For instance for temperature, there is a cold, normal and warm class. Commonly multiple ensemble forecasting systems are available but some models may be more competitive than others for the variable, region and user need under consideration. Therefore, when calculating the category probabilities, the ensemble members of the different forecasting system may be differently weighted. Such weighting is typically done by comparison of the ensemble forecasts with observations.

Description of the tool: The tool considers all forecasts (all members from all forecasting systems) and converts them into occurrence probabilities of different categories. The amount of categories can be changed and are taken as the climatological quantiles (e.g. terciles), extracted from the observational data. The methods that are available to combine the ensemble forecasting models into probabilistic categorical forecasts are: 1) ensemble pooling where all ensemble members of all ensemble systems are weighted equally, 2) model combination where each model system is weighted equally, and, 3) model weighting. The model weighting method is described in Rajagopalan et al. (2002), Robertson et al. 2004 and Van Schaeybroeck and Vannitsem (2019). More specifically, this method uses different weights for the occurence probability predicted by the available models and by a climatological model and optimizes the weights by minimizing the ignorance score. Finally, the function can also be used to categorize the observations in the categorical quantiles.

Usage

CST_CategoricalEnsCombination(
  exp,
  obs,
  cat.method = "pool",
  eval.method = "leave-k-out",
  amt.cat = 3,
  k = 1,
  tail.out = TRUE,
  ...
)

Arguments

Value

An object of class s2dv_cube containing the categorical forecasts in the element called $data. The first two dimensions of the returned object are named dataset and member and are both of size one. An additional dimension named category is introduced and is of size amt.cat.

Author(s)

Bert Van Schaeybroeck, bertvs@meteo.be

References

Rajagopalan, B., Lall, U., & Zebiak, S. E. (2002). Categorical climate forecasts through regularization and optimal combination of multiple GCM ensembles. Monthly Weather Review, 130(7), 1792-1811.

Robertson, A. W., Lall, U., Zebiak, S. E., & Goddard, L. (2004). Improved combination of multiple atmospheric GCM ensembles for seasonal prediction. Monthly Weather Review, 132(12), 2732-2744.

Van Schaeybroeck, B., & Vannitsem, S. (2019). Postprocessing of Long-Range Forecasts. In Statistical Postprocessing of Ensemble Forecasts (pp. 267-290).

Examples

mod1 <- 1 : (2 * 2* 4 * 5 * 2 * 2)
dim(mod1) <- c(dataset = 2, member = 2, sdate = 4, ftime = 5, lat = 2, lon = 2)
mod1[2, 1, , , , ] <- NA
datasets <- c("MF", "UKMO")
obs1 <- 1 : (1 * 1 * 4 * 5 * 2 * 2)
dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 2, lon = 2)
lon <- seq(0, 30, 5)
lat <- seq(0, 25, 5)
coords <- list(lat = lat, lon = lon)
attrs <- list(Datasets = datasets)
exp <- list(data = mod1, coords = coords, attrs = attrs)
obs <- list(data = obs1, coords = coords)
attr(exp, 'class') <- 's2dv_cube'
attr(obs, 'class') <- 's2dv_cube'
a <- CST_CategoricalEnsCombination(exp = exp, obs = obs, amt.cat = 3, 
                                  cat.method = "mmw") 

Change the name of one or more dimensions for an object of class s2dv_cube

Description

Change the names of the dimensions specified in 'original_names' to the names in 'new_names'. The coordinate names and the dimensions of any attributes are also modified accordingly.

Usage

CST_ChangeDimNames(data, original_names, new_names)

Arguments

Value

An object of class s2dv_cube with similar data, coordinates and attributes as the data input, but with modified dimension names.

Author(s)

Agudetse Roures Victoria, victoria.agudetse@bsc.es

Examples

# Example with sample data:
# Check original dimensions and coordinates
lonlat_temp$exp$dims
names(lonlat_temp$exp$coords)
dim(lonlat_temp$exp$attrs$Dates)
# Change 'dataset' to 'dat' and 'ftime' to 'time'
exp <- CST_ChangeDimNames(lonlat_temp$exp,
                         original_names = c("dataset", "ftime"),
                         new_names = c("dat", "time"))
# Check new dimensions and coordinates
exp$dims
names(exp$coords)
dim(exp$attrs$Dates)

Performing a Bias Correction conditioned by the dynamical properties of the data.

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'

Usage

CST_DynBiasCorrection(
  exp,
  obs,
  method = "QUANT",
  wetday = FALSE,
  proxy = "dim",
  quanti,
  ncores = NULL
)

Arguments

Value

dynbias An s2dvcube object with a bias correction performed conditioned by local dimension 'dim' or inverse of persistence 'theta'.

Author(s)

Carmen Alvarez-Castro, carmen.alvarez-castro@cmcc.it

Maria M. Chaves-Montero, mdm.chaves-montero@cmcc.it

Veronica Torralba, veronica.torralba@cmcc.it

Davide Faranda, davide.faranda@lsce.ipsl.fr

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: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.

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)
time_obsL <- as.POSIXct(paste(rep("01", 100), rep("01", 100), 1920:2019, sep = "-"), 
                       format = "%d-%m-%y")
time_expL <- as.POSIXct(paste(rep("01", 100), rep("01", 100), 1929:2019, sep = "-"), 
                       format = "%d-%m-%y")
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, coords = list(lon = lon, lat = lat),
                 Dates = time_expL)
obsL <- s2dv_cube(data = obsL, coords = list(lon = lon, lat = lat),
                 Dates = 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)

Ensemble clustering

Description

This function performs a clustering on members/starting dates and returns a number of scenarios, with representative members for each of them. The clustering is performed in a reduced EOF space.

Motivation: Ensemble forecasts give a probabilistic insight of average weather conditions on extended timescales, i.e. from sub-seasonal to seasonal and beyond. With large ensembles, it is often an advantage to be able to group members according to similar characteristics and to select the most representative member for each cluster. This can be useful to characterize the most probable forecast scenarios in a multi-model (or single model) ensemble prediction. This approach, applied at a regional level, can also be used to identify the subset of ensemble members that best represent the full range of possible solutions for downscaling applications. The choice of the ensemble members is made flexible in order to meet the requirements of specific (regional) climate information products, to be tailored for different regions and user needs.

Description of the tool: EnsClustering is a cluster analysis tool, based on the k-means algorithm, for ensemble predictions. The aim is to group ensemble members according to similar characteristics and to select the most representative member for each cluster. The user chooses which feature of the data is used to group the ensemble members by clustering: time mean, maximum, a certain percentile (e.g., 75 time period. For each ensemble member this value is computed at each grid point, obtaining N lat-lon maps, where N is the number of ensemble members. The anomaly is computed subtracting the ensemble mean of these maps to each of the single maps. The anomaly is therefore computed with respect to the ensemble members (and not with respect to the time) and the Empirical Orthogonal Function (EOF) analysis is applied to these anomaly maps. Regarding the EOF analysis, the user can choose either how many Principal Components (PCs) to retain or the percentage of explained variance to keep. After reducing dimensionality via EOF analysis, k-means analysis is applied using the desired subset of PCs.

The major final outputs are the classification in clusters, i.e. which member belongs to which cluster (in k-means analysis the number k of clusters needs to be defined prior to the analysis) and the most representative member for each cluster, which is the closest member to the cluster centroid. Other outputs refer to the statistics of clustering: in the PC space, the minimum and the maximum distance between a member in a cluster and the cluster centroid (i.e. the closest and the furthest member), the intra-cluster standard deviation for each cluster (i.e. how much the cluster is compact).

Usage

CST_EnsClustering(
  exp,
  time_moment = "mean",
  numclus = NULL,
  lon_lim = NULL,
  lat_lim = NULL,
  variance_explained = 80,
  numpcs = NULL,
  time_dim = NULL,
  time_percentile = 90,
  cluster_dim = "member",
  verbose = FALSE
)

Arguments

Value

A list with elements $cluster (cluster assigned for each member), $freq (relative frequency of each cluster), $closest_member (representative member for each cluster), $repr_field (list of fields for each representative member), composites (list of mean fields for each cluster), $lon (selected longitudes of output fields), $lat (selected longitudes of output fields).

Author(s)

Federico Fabiano - ISAC-CNR, f.fabiano@isac.cnr.it

Ignazio Giuntoli - ISAC-CNR, i.giuntoli@isac.cnr.it

Danila Volpi - ISAC-CNR, d.volpi@isac.cnr.it

Paolo Davini - ISAC-CNR, p.davini@isac.cnr.it

Jost von Hardenberg - ISAC-CNR, j.vonhardenberg@isac.cnr.it

Examples

dat_exp <- array(abs(rnorm(1152))*275, dim = c(dataset = 1, member = 4, 
                                              sdate = 6, ftime = 3, 
                                              lat = 4, lon = 4))
lon <- seq(0, 3)
lat <- seq(48, 45)
coords <- list(lon = lon, lat = lat)
exp <- list(data = dat_exp, coords = coords)
attr(exp, 'class') <- 's2dv_cube'
res <- CST_EnsClustering(exp = exp, numclus = 3,
                        cluster_dim = c("sdate"))

Add a named dimension to an object of class s2dv_cube

Description

Insert an extra dimension into an array at position 'posdim' with length 'lendim'. The array in data repeats along the new dimension. The dimensions, coordinates and attributes are modified accordingly.

Usage

CST_InsertDim(data, posdim, lendim, name, values = NULL)

Arguments

Value

An object of class s2dv_cube with similar data, coordinates and attributes as the data input, but with an additional dimension.

Author(s)

Agudetse Roures Victoria, victoria.agudetse@bsc.es

See Also

InsertDim

Examples

#Example with sample data:
# Check original dimensions and coordinates
lonlat_temp$exp$dims
names(lonlat_temp$exp$coords)
# Add 'variable' dimension
exp <- CST_InsertDim(lonlat_temp$exp,
                    posdim = 2,
                    lendim = 1,
                    name = "variable",
                    values = c("tas"))
# Check new dimensions and coordinates
exp$dims
exp$coords$variable

CSTools Data Retreival Function

Description

This function aggregates, subsets and retrieves sub-seasonal, seasonal, decadal or climate projection data from NetCDF files in a local file system or on remote OPeNDAP servers, and arranges it for easy application of the CSTools functions.

Usage

CST_Load(...)

Arguments

Details

It receives any number of parameters ('...') that are automatically forwarded to the 's2dv::Load' function. See details in '?s2dv::Load'.

It is recommended to use this function in combination with the 'zeallot::"

Value

A list with one or two S3 objects, named 'exp' and 'obs', of the class 's2dv_cube', containing experimental and date-corresponding observational data, respectively. These 's2dv_cube's can be ingested by other functions in CSTools. If the parameter ‘exp' in the call to 'CST_Load' is set to 'NULL', then only the ’obs' component is returned, and viceversa.

Author(s)

Nicolau Manubens, nicolau.manubens@bsc.es

Examples

## Not run: 
library(zeallot)
startDates <- c('20001101', '20011101', '20021101', 
                '20031101', '20041101', '20051101')
c(exp, obs) %<-% 
  CST_Load(
    var = 'tas', 
    exp = 'system5c3s', 
    obs = 'era5', 
    nmember = 15,
    sdates = startDates,
    leadtimemax = 3,
    latmin = 27, latmax = 48,
    lonmin = -12, lonmax = 40,
    output = 'lonlat',
    nprocs = 1
  )

## End(Not run)

Function to Merge Dimensions

Description

This function merges two dimensions of the array data in a 's2dv_cube' object into one. The user can select the dimensions to merge and provide the final name of the dimension. The user can select to remove NA values or keep them.

Usage

CST_MergeDims(
  data,
  merge_dims = c("ftime", "monthly"),
  rename_dim = NULL,
  na.rm = FALSE
)

Arguments

Value

An object of class 's2dv_cube' with the specified dimensions merged into a single dimension.

Author(s)

Nuria Perez-Zanon, nuria.perez@bsc.es

Examples

data <- 1 : c(2 * 3 * 4 * 5 * 6 * 7)
dim(data) <- c(time = 7, lat = 2, lon = 3, monthly = 4, member = 6,
              dataset = 5, var = 1)
data[2,,,,,,] <- NA
data[c(3,27)] <- NA
data <- list(data = data)
class(data) <- 's2dv_cube'
new_data <- CST_MergeDims(data, merge_dims = c('time', 'monthly'))
new_data <- CST_MergeDims(data, merge_dims = c('lon', 'lat'), rename_dim = 'grid')
new_data <- CST_MergeDims(data, merge_dims = c('time', 'monthly'), na.rm = TRUE)

EOF analysis of multiple variables

Description

This function performs EOF analysis over multiple variables, accepting in input a list of CSTools objects. Based on Singular Value Decomposition. For each field the EOFs are computed and the corresponding PCs are standardized (unit variance, zero mean); the minimum number of principal components needed to reach the user-defined variance is retained. The function weights the input data for the latitude cosine square root.

Usage

CST_MultiEOF(
  datalist,
  lon_dim = "lon",
  lat_dim = "lat",
  time_dim = "ftime",
  sdate_dim = "sdate",
  var_dim = "var",
  neof_max = 40,
  neof_composed = 5,
  minvar = 0.6,
  lon_lim = NULL,
  lat_lim = NULL,
  ncores = NULL
)

Arguments

Value

A list containing:

Author(s)

Jost von Hardenberg - ISAC-CNR, j.vonhardenberg@isac.cnr.it

Paolo Davini - ISAC-CNR, p.davini@isac.cnr.it

Examples

seq <- 1 : (2 * 3 * 4 * 5 * 6 * 8)
mod1 <- sin( 0.7 + seq )^2 + cos( seq ^ 2 * 1.22  )
dim(mod1) <- c(dataset = 2, member = 3, sdate = 4, ftime = 5, lat = 6, 
              lon = 8)
mod2 <- sin( seq * 2 ) ^ 3 + cos( seq ^ 2 )
dim(mod2) <- c(dataset = 2, member = 3, sdate = 4, ftime = 5, lat = 6, 
              lon = 8)
lon <- seq(0, 35, 5)
lat <- seq(0, 25, 5)
exp1 <- list(data = mod1, coords = list(lat = lat, lon = lon))
exp2 <- list(data = mod2, coords = list(lat = lat, lon = lon))
attr(exp1, 'class') <- 's2dv_cube'
attr(exp2, 'class') <- 's2dv_cube'
d = as.POSIXct(c("2017/01/01", "2017/01/02", "2017/01/03", "2017/01/04", 
                "2017/01/05", "2018/01/01", "2018/01/02", "2018/01/03",
                "2018/01/04", "2018/01/05", "2019/01/01", "2019/01/02", 
                "2019/01/03", "2019/01/04", "2019/01/05", "2020/01/01", 
                "2020/01/02", "2020/01/03", "2020/01/04", "2020/01/05"))
exp1$attrs$Dates = d
exp2$attrs$Dates = d

cal <- CST_MultiEOF(datalist = list(exp1, exp2), neof_composed = 2)

Multiple Metrics applied in Multiple Model Anomalies

Description

This function calculates correlation (Anomaly Correlation Coefficient; ACC), root mean square error (RMS) and the root mean square error skill score (RMSSS) of individual anomaly models and multi-models mean (if desired) with the observations.

Usage

CST_MultiMetric(
  exp,
  obs,
  metric = "correlation",
  multimodel = TRUE,
  time_dim = "ftime",
  memb_dim = "member",
  sdate_dim = "sdate"
)

Arguments

Value

An object of class s2dv_cube containing the statistics of the selected metric in the element $data which is a list of arrays: for the metric requested and others for statistics about its signeificance. The arrays have two dataset dimensions equal to the 'dataset' dimension in the exp$data and obs$data inputs. If multimodel is TRUE, the first position in the first 'nexp' dimension correspons to the Multi-Model Mean.

Author(s)

Mishra Niti, niti.mishra@bsc.es

Perez-Zanon Nuria, nuria.perez@bsc.es

References

Mishra, N., Prodhomme, C., & Guemas, V. (n.d.). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29-31. doi:10.1007/s00382-018-4404-z

See Also

Corr, RMS, RMSSS and CST_Load

Examples

mod <- rnorm(2*2*4*5*2*2)
dim(mod) <- c(dataset = 2, member = 2, sdate = 4, ftime = 5, lat = 2, lon = 2)
obs <- rnorm(1*1*4*5*2*2)
dim(obs) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 2, lon = 2)
lon <- seq(0, 30, 5)
lat <- seq(0, 25, 5)
coords <- list(lat = lat, lon = lon)
exp <- list(data = mod, coords = coords)
obs <- list(data = obs, coords = coords)
attr(exp, 'class') <- 's2dv_cube'
attr(obs, 'class') <- 's2dv_cube'
a <- CST_MultiMetric(exp = exp, obs = obs)

Multivariate Root Mean Square Error (RMSE)

Description

This function calculates the RMSE from multiple variables, as the mean of each variable's RMSE scaled by its observed standard deviation. Variables can be weighted based on their relative importance (defined by the user).

Usage

CST_MultivarRMSE(
  exp,
  obs,
  weight = NULL,
  memb_dim = "member",
  dat_dim = "dataset",
  sdate_dim = "sdate",
  ftime_dim = "ftime"
)

Arguments

Value

An object of class s2dv_cube containing the RMSE in the element $data which is an array with two datset dimensions equal to the 'dataset' dimension in the exp$data and obs$data inputs. An array with dimensions: c(number of exp, number of obs, 1 (the multivariate RMSE value), number of lat, number of lon)

Author(s)

Deborah Verfaillie, deborah.verfaillie@bsc.es

See Also

RMS and CST_Load

Examples

# Example with 2 variables
mod1 <- abs(rnorm(1 * 3 * 4 * 5 * 6 * 7))
mod2 <- abs(rnorm(1 * 3 * 4 * 5 * 6 * 7))
dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
dim(mod2) <- c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
obs1 <- abs(rnorm(1 * 1 * 4 * 5 * 6 * 7))
obs2 <- abs(rnorm(1 * 1 * 4 * 5 * 6 * 7))
dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
dim(obs2) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
lon <- seq(0, 30, 5)
lat <- seq(0, 25, 5)
coords <- list(lat = lat, lon = lon)
exp1 <- list(data = mod1, coords = coords, 
            attrs = list(Datasets = "EXP1", source_files = "file1", 
                         Variable = list(varName = 'pre')))
exp2 <- list(data = mod2, coords = coords, 
            attrs = list(Datasets = "EXP2", source_files = "file2", 
                         Variable = list(varName = 'tas')))
obs1 <- list(data = obs1, coords = coords, 
            attrs = list(Datasets = "OBS1", source_files = "file1", 
                         Variable = list(varName = 'pre')))
obs2 <- list(data = obs2, coords = coords, 
            attrs = list(Datasets = "OBS2", source_files = "file2", 
                         Variable = list(varName = 'tas')))
attr(exp1, 'class') <- 's2dv_cube'
attr(exp2, 'class') <- 's2dv_cube'
attr(obs1, 'class') <- 's2dv_cube'
attr(obs2, 'class') <- 's2dv_cube'
anom1 <- CST_Anomaly(exp1, obs1, cross = TRUE, memb = TRUE)
anom2 <- CST_Anomaly(exp2, obs2, cross = TRUE, memb = TRUE)
ano_exp <- list(anom1$exp, anom2$exp)
ano_obs <- list(anom1$obs, anom2$obs)
a <- CST_MultivarRMSE(exp = ano_exp, obs = ano_obs, weight = c(1, 2))

Computing two dinamical proxies of the attractor in s2dv_cube.

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 Function based on the matlab code (davide.faranda@lsce.ipsl.fr) used in

Usage

CST_ProxiesAttractor(data, quanti, ncores = NULL)

Arguments

Value

dim and theta

Author(s)

Carmen Alvarez-Castro, carmen.alvarez-castro@cmcc.it

Maria M. Chaves-Montero, mdm.chaves-montero@cmcc.it

Veronica Torralba, veronica.torralba@cmcc.it

Davide Faranda, davide.faranda@lsce.ipsl.fr

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: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.

Examples

# Example 1: Computing the attractor using simple s2dv data
obs <- rnorm(2 * 3 * 4 * 8 * 8)
dim(obs) <- c(dataset = 1, member = 2, sdate = 3, ftime = 4, lat = 8, lon = 8)
lon <- seq(10, 13.5, 0.5)
lat <- seq(40, 43.5, 0.5)
coords <- list(lon = lon, lat = lat)
data <- list(data = obs, coords = coords)
class(data) <- "s2dv_cube"
attractor <- CST_ProxiesAttractor(data = data, quanti = 0.6)

Quantile Mapping for seasonal or decadal forecast data

Description

This function is a wrapper of fitQmap and doQmap from package 'qmap' to be applied on the object of class 's2dv_cube'. The quantile mapping adjustment between an experiment, typically a hindcast, and observation is applied to the experiment itself or to a provided forecast.

Usage

CST_QuantileMapping(
  exp,
  obs,
  exp_cor = NULL,
  sdate_dim = "sdate",
  memb_dim = "member",
  window_dim = NULL,
  method = "QUANT",
  na.rm = FALSE,
  eval.method = "leave-k-out",
  k = 1,
  ncores = NULL,
  ...
)

Arguments

Value

An object of class s2dv_cube containing the experimental data after applying the quantile mapping correction.

Author(s)

Nuria Perez-Zanon, nuria.perez@bsc.es

See Also

fitQmap and doQmap

Examples

# Use synthetic data
exp <- NULL
exp$data <- 1 : c(1 * 3 * 5 * 4 * 3 * 2)
dim(exp$data) <- c(dataset = 1, member = 3, sdate = 5, ftime = 4,
                  lat = 3, lon = 2)
class(exp) <- 's2dv_cube'
obs <- NULL
obs$data <- 101 : c(100 + 1 * 1 * 5 * 4 * 3 * 2)
dim(obs$data) <- c(dataset = 1, member = 1, sdate = 5, ftime = 4,
                  lat = 3, lon = 2)
class(obs) <- 's2dv_cube'
res <- CST_QuantileMapping(exp, obs)

RainFARM spectral slopes from a CSTools object

Description

This function computes spatial spectral slopes from a CSTools object to be used for RainFARM stochastic precipitation downscaling method and accepts a CSTools object (of the class 's2dv_cube') as input.

Usage

CST_RFSlope(data, kmin = 1, time_dim = NULL, ncores = NULL)

Arguments

Value

CST_RFSlope() returns spectral slopes using the RainFARM convention (the logarithmic slope of k*|A(k)|^2 where A(k) are the spectral amplitudes). The returned array has the same dimensions as the exp element of the input object, minus the dimensions specified by lon_dim, lat_dim and time_dim.

Author(s)

Jost von Hardenberg - ISAC-CNR, j.vonhardenberg@isac.cnr.it

Examples

exp <- 1 : (2 * 3 * 4 * 8 * 8)
dim(exp) <- c(dataset = 1, member = 2, sdate = 3, ftime = 4, lat = 8, lon = 8)
lon <- seq(10, 13.5, 0.5)
lat <- seq(40, 43.5, 0.5)
coords <- list(lon = lon, lat = lat)
data <- list(data = exp, coords = coords)
class(data) <- 's2dv_cube'
slopes <- CST_RFSlope(data)

Temperature downscaling of a CSTools object using lapse rate correction or a reference field

Description

This function implements a simple lapse rate correction of a temperature field (an object of class 's2dv_cube' as provided by 'CST_Load') as input. The input lon grid must be increasing (but can be modulo 360). The input lat grid can be irregularly spaced (e.g. a Gaussian grid) The output grid can be irregularly spaced in lon and/or lat.

Usage

CST_RFTemp(
  data,
  oro,
  xlim = NULL,
  ylim = NULL,
  lapse = 6.5,
  lon_dim = "lon",
  lat_dim = "lat",
  time_dim = NULL,
  nolapse = FALSE,
  verbose = FALSE,
  compute_delta = FALSE,
  method = "bilinear",
  delta = NULL
)

Arguments

Value

CST_RFTemp() returns a downscaled CSTools object (i.e., of the class 's2dv_cube').

Author(s)

Jost von Hardenberg - ISAC-CNR, j.vonhardenberg@isac.cnr.it

References

Method described in ERA4CS MEDSCOPE milestone M3.2: High-quality climate prediction data available to WP4 here: https://www.medscope-project.eu/the-project/deliverables-reports/ and in H2020 ECOPOTENTIAL Deliverable No. 8.1: High resolution (1-10 km) climate, land use and ocean change scenarios available here: https://ec.europa.eu/research/participants/documents/downloadPublic?documentIds=080166e5b6cd2324&appId=PPGMS

Examples

# Generate simple synthetic data and downscale by factor 4
t <- rnorm(7 * 6 * 2 * 3 * 4)*10 + 273.15 + 10
dim(t) <- c(dataset = 1, member = 2, sdate = 3, ftime = 4, lat = 6, lon = 7)
lon <- seq(3, 9, 1)
lat <- seq(42, 47, 1)
coords <- list(lat = lat, lon = lon)
exp <- list(data = t, coords = coords)
attr(exp, 'class') <- 's2dv_cube'
o <- runif(29*29)*3000
dim(o) <- c(lats = 29, lons = 29)
lon <- seq(3, 10, 0.25)
lat <- seq(41, 48, 0.25)
coords <- list(lat = lat, lon = lon)
oro <- list(data = o, coords = coords)
attr(oro, 'class') <- 's2dv_cube'
res <- CST_RFTemp(data = exp, oro = oro, xlim = c(4,8), ylim = c(43, 46), 
                 lapse = 6.5, time_dim = 'ftime',
                 lon_dim = 'lon', lat_dim = 'lat')

Compute climatological weights for RainFARM stochastic precipitation downscaling

Description

Compute climatological ("orographic") weights from a fine-scale precipitation climatology file.

Usage

CST_RFWeights(
  climfile,
  nf,
  lon,
  lat,
  varname = NULL,
  fsmooth = TRUE,
  lonname = "lon",
  latname = "lat",
  ncores = NULL
)

Arguments

Value

An object of class 's2dv_cube' containing in matrix data the weights with dimensions (lon, lat).

Author(s)

Jost von Hardenberg - ISAC-CNR, j.vonhardenberg@isac.cnr.it

References

Terzago, S., Palazzi, E., & von Hardenberg, J. (2018). Stochastic downscaling of precipitation in complex orography: A simple method to reproduce a realistic fine-scale climatology. Natural Hazards and Earth System Sciences, 18(11), 2825-2840. doi:10.5194/nhess-18-2825-2018.

Examples

# Create weights to be used with the CST_RainFARM() or RainFARM() functions
# using an external random data in the form of 's2dv_cube'.
obs <- rnorm(2 * 3 * 4 * 8 * 8)
dim(obs) <- c(dataset = 1, member = 2, sdate = 3, ftime = 4, lat = 8, lon = 8)
lon <- seq(10, 13.5, 0.5)
lat <- seq(40, 43.5, 0.5)
coords <- list(lon = lon, lat = lat)
data <- list(data = obs, coords = coords)
class(data) <- "s2dv_cube"
res <- CST_RFWeights(climfile = data, nf = 3, lon, lat, lonname = 'lon', 
                    latname = 'lat', fsmooth = TRUE)

RainFARM stochastic precipitation downscaling of a CSTools object

Description

This function implements the RainFARM stochastic precipitation downscaling method and accepts a CSTools object (an object of the class 's2dv_cube' as provided by 'CST_Load') as input. Adapted for climate downscaling and including orographic correction as described in Terzago et al. 2018.

Usage

CST_RainFARM(
  data,
  weights = 1,
  slope = 0,
  nf,
  kmin = 1,
  nens = 1,
  fglob = FALSE,
  fsmooth = TRUE,
  nprocs = 1,
  time_dim = NULL,
  verbose = FALSE,
  drop_realization_dim = FALSE
)

Arguments

Details

Wether parameter 'slope' and 'weights' presents seasonality dependency, a dimension name should match between these parameters and the input data in parameter 'data'. See example 2 below where weights and slope vary with 'sdate' dimension.

Value

CST_RainFARM() returns a downscaled CSTools object (i.e., of the class 's2dv_cube'). If nens > 1 an additional dimension named "realization" is added to the $data array after the "member" dimension (unless drop_realization_dim = TRUE is specified). The ordering of the remaining dimensions in the $data element of the input object is maintained.

Author(s)

Jost von Hardenberg - ISAC-CNR, j.vonhardenberg@isac.cnr.it

References

Terzago, S. et al. (2018). NHESS 18(11), 2825-2840. doi:10.5194/nhess-18-2825-2018; D'Onofrio et al. (2014), J of Hydrometeorology 15, 830-843; Rebora et. al. (2006), JHM 7, 724.

Examples

# Example 1: using CST_RainFARM for a CSTools object
nf <- 8   # Choose a downscaling by factor 8
exp <- 1 : (2 * 3 * 4 * 8 * 8)
dim(exp) <- c(dataset = 1, member = 2, sdate = 3, ftime = 4, lat = 8, lon = 8)
lon <- seq(10, 13.5, 0.5)
lat <- seq(40, 43.5, 0.5)
coords <- list(lon = lon, lat = lat)
data <- list(data = exp, coords = coords)
class(data) <- 's2dv_cube'
# Create a test array of weights
ww <- array(1., dim = c(lon = 8 * nf, lat = 8 * nf))
res <- CST_RainFARM(data, nf = nf, weights = ww, nens = 3, time_dim = 'ftime')

Function for matching a field of anomalies with a set of maps used as a reference (e.g. clusters obtained from the WeatherRegime function)

Description

This function performs the matching between a field of anomalies and a set of maps which will be used as a reference. The anomalies will be assigned to the reference map for which the minimum Eucledian distance (method =’distance’) or highest spatial correlation (method = 'ACC') is obtained.

Usage

CST_RegimesAssign(
  data,
  ref_maps,
  method = "distance",
  composite = FALSE,
  memb = FALSE,
  ncores = NULL
)

Arguments

Value

A list with two elements $data (a 's2dv_cube' object containing the composites cluster=1,..,K for case (*1) or only k=1 for any specific cluster, i.e., case (*2)) (only when composite = 'TRUE') and $statistics that includes $pvalue (array with the same structure as $data containing the pvalue of the composites obtained through a t-test that accounts for the serial dependence of the data with the same structure as Composite.)(only when composite = 'TRUE'), $cluster (array with the same dimensions as data (except latitude and longitude which are removed) indicating the ref_maps to which each point is allocated.), $frequency (A vector of integers (from k=1,...k n reference maps) indicating the percentage of assignations corresponding to each map.).

Author(s)

Verónica Torralba - BSC, veronica.torralba@bsc.es

References

Torralba, V. (2019) Seasonal climate prediction for the wind energy sector: methods and tools for the development of a climate service. Thesis. Available online: https://eprints.ucm.es/56841/

Examples

data <- array(abs(rnorm(1280, 282.7, 6.4)), dim = c(dataset = 2, member = 2, 
                                                   sdate = 3, ftime = 3, 
                                                   lat = 4, lon = 4))
coords <- list(lon = seq(0, 3), lat = seq(47, 44))
exp <- list(data = data, coords = coords)
class(exp) <- 's2dv_cube'
regimes <- CST_WeatherRegimes(data = exp, EOFs = FALSE, 
                             ncenters = 4)
res1 <- CST_RegimesAssign(data = exp, ref_maps = regimes, 
                         composite = FALSE)

Save objects of class 's2dv_cube' to data in NetCDF format

Description

This function allows to divide and save a object of class 's2dv_cube' into a NetCDF file, allowing to reload the saved data using CST_Start or CST_Load functions. It also allows to save any 's2dv_cube' object that follows the NetCDF attributes conventions.

Usage

CST_SaveExp(
  data,
  destination = tempdir(),
  startdates = NULL,
  sdate_dim = "sdate",
  ftime_dim = "time",
  memb_dim = "member",
  dat_dim = "dataset",
  var_dim = "var",
  drop_dims = NULL,
  single_file = FALSE,
  extra_string = NULL,
  global_attrs = NULL,
  units_hours_since = FALSE
)

Arguments

Value

Multiple or single NetCDF files containing the data array.

Author(s)

Perez-Zanon Nuria, nuria.perez@bsc.es

See Also

Start, as.s2dv_cube and s2dv_cube

Examples

data <- lonlat_temp_st$exp
CST_SaveExp(data = data, ftime_dim = 'ftime', var_dim = 'var', 
           dat_dim = 'dataset', sdate_dim = 'sdate')

Function to Split Dimension

Description

This function split a dimension in two. The user can select the dimension to split and provide indices indicating how to split that dimension or dates and the frequency expected (monthly or by day, month and year). The user can also provide a numeric frequency indicating the length of each division.

Usage

CST_SplitDim(
  data,
  split_dim = "time",
  indices = NULL,
  freq = "monthly",
  new_dim_name = NULL,
  insert_ftime = NULL,
  ftime_dim = "time",
  sdate_dim = "sdate",
  return_indices = FALSE
)

Arguments

Details

Parameter 'insert_ftime' has been included for the case of using daily data, requiring split the temporal dimensions by months (or similar) and the first lead time doesn't correspondt to the 1st day of the month. In this case, the insert_ftime could be used, to get a final output correctly organized. E.g.: leadtime 1 is the 2nd of November and the input time series extend to the 31st of December. When requiring split by month with inset_ftime = 1, the 'monthly' dimension of length two will indicate the month (position 1 for November and position 2 for December), dimension 'time' will be length 31. For November, the position 1 and 31 will be NAs, while from positon 2 to 30 will be filled with the data provided. This allows to select correctly days trhough time dimension.

Value

An object of class 's2dv_cube' with the specified dimension split into a new dimension, preserving all other original dimensions, coordinates, and attributes (including Dates).

Author(s)

Nuria Perez-Zanon, nuria.perez@bsc.es

Examples

data <- 1 : 20
dim(data) <- c(time = 10, lat = 2)
data <-list(data = data)
class(data) <- 's2dv_cube'
indices <- c(rep(1,5), rep(2,5))
new_data <- CST_SplitDim(data, indices = indices)
time <- c(seq(ISOdate(1903, 1, 1), ISOdate(1903, 1, 4), "days"),
         seq(ISOdate(1903, 2, 1), ISOdate(1903, 2, 4), "days"),
         seq(ISOdate(1904, 1, 1), ISOdate(1904, 1, 2), "days"))
data <- list(data = data$data, Dates = time)
class(data) <- 's2dv_cube'
new_data <- CST_SplitDim(data, indices = time)
new_data <- CST_SplitDim(data, indices = time, freq = 'day')
new_data <- CST_SplitDim(data, indices = time, freq = 'month')
new_data <- CST_SplitDim(data, indices = time, freq = 'year')

CSTools data retrieval function using Start

Description

This function aggregates, subsets and retrieves sub-seasonal, seasonal, decadal or climate projection data from NetCDF files in a local file system and arranges it for easy application of the CSTools functions. It calls the function Start from startR, which is an R package started at BSC with the aim to develop a tool that allows the user to automatically process large multidimensional distributed data sets. Then, the output is transformed into 's2dv_cube' object.

Usage

CST_Start(...)

Arguments

Details

It receives any number of parameters ('...') that are automatically forwarded to the 'startR::Start' function. See details in '?startR::Start'. The auxiliary startR functions (e.g. indices(), values(), Sort(), CircularSort(), CDORemapper(), ...) can be used to define dimensions.

CST_Start() uses as.s2dv_cube() to transform the output into an s2dv_cube object. The as.s2dv_cube() function is designed to be used with data that has been retrieved into memory. To avoid errors, please ensure that CST_Start(..., retrieve = TRUE) is specified.

Value

An object of class 's2dv_cube' containing the subsetted and aggregated data retrieved from the NetCDF files using startR. The data structure is compatible with CSTools functions.

Examples

## Not run: 
 sdates <- c('20101101', '20111101', '20121101')
 latmin <- 44
 latmax <- 47
 lonmin <- 6
 lonmax <- 9
 data <- CST_Start(dat = path,
                   var = 'prlr',
                   ensemble = indices(1:6),
                   sdate = sdates,
                   time = 121:151,
                   latitude = values(list(latmin, latmax)),
                   longitude = values(list(lonmin, lonmax)),
                   synonims = list(longitude = c('lon', 'longitude'),
                                   latitude = c('lat', 'latitude')),
                   return_vars = list(time = 'sdate',
                                      longitude = NULL, latitude = NULL),
                   retrieve = TRUE)

## End(Not run) 

Subset an object of class s2dv_cube

Description

This function allows to subset (i.e. slice, take a chunk of) the data inside an object of class s2dv_cube and modify the dimensions, coordinates and attributes accordingly, removing any variables, time steps and spatial coordinates that are dropped when subsetting. It ensures that the information inside the s2dv_cube remains coherent with the data it contains.

As in the function Subset from the ClimProjDiags package, the dimensions to subset along can be specified via the parameter along either with integer indices or by their name.

There are additional ways to adjust which dimensions are dropped in the resulting object: either to drop all, to drop none, to drop only the ones that have been sliced or to drop only the ones that have not been sliced.

The load_parameters and when attributes of the original cube are preserved. The source_files attribute is subset along the var_dim and dat_dim dimensions.

Usage

CST_Subset(x, along, indices, drop = FALSE, var_dim = NULL, dat_dim = NULL)

Arguments

Value

An object of class s2dv_cube with similar data, coordinates and attributes as the x input, but with trimmed or dropped dimensions.

Author(s)

Agudetse Roures Victoria, victoria.agudetse@bsc.es

See Also

Subset

Examples

#Example with sample data:
# Check original dimensions and coordinates
lonlat_temp$exp$dims
names(lonlat_temp$exp$coords)
# Subset the s2dv_cube
exp_subset <- CST_Subset(lonlat_temp$exp,
                        along = c("lat", "lon"),
                        indices = list(1:10, 1:10),
                        drop = 'non-selected')
# Check new dimensions and coordinates
exp_subset$dims
names(exp_subset$coords)

Generate a Summary of the data and metadata in the s2dv_cube object

Description

This function prints the summary of the data and metadata of an object of class s2dv_cube.

Usage

CST_Summary(data, show_loaded_files = FALSE, show_NA = FALSE, var_dim = "var")

Arguments

Details

The function uses the data and metadata from the loaded s2dv_cube object to generate a summary of the object.The summary includes : - months: Months that have been loaded. - range: Range of the dates that have been loaded. - dimensions: Object dimensions. - Statistical summary of the data in data: Basic statistical summary of the data. - Variable: Loaded Variables, along with their units (units:) - NA-Indices per Dimension: Index with NA values per dimension - Number of NAs in NA-Indices per Dimensions: Number of NAs, in the Indices with NA values per dimension - Loaded files: Successfully loaded Files

Value

Does not return a value. This function prints a detailed summary of an s2dv_cube object to the console.

Author(s)

Theertha Kariyathan, theertha.kariyathan@bsc.es

See Also

CST_Start or s2dv_cube for creating an s2dv cube object.

Examples

# Example 1:
CST_Summary(data = lonlat_temp_st$exp)

# Example 2:
## Not run: 
# s2dv cube paths
repos1 <- "/esarchive/exp/ecmwf/system4_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc"
repos2 <- "/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc"

# Create data cube
data <- CST_Start(dat = list(
                  list(name = 'system4_m1', path = repos1),
                  list(name = 'system5_m1', path = repos2)),
                  var = c('tas', 'sfcWind'),
                  sdate = '20170101',
                  ensemble = indices(1),
                  time = indices(1:3),
                  lat = indices(1:5),
                  lon = indices(1:5),
                  synonims = list(lat = c('lat', 'latitude'),
                                  lon = c('lon', 'longitude')),
                  return_vars = list(time = 'sdate',
                                     longitude = 'dat',
                                     latitude = 'dat'),
                  metadata_dims = c('dat', 'var'),
                  retrieve = TRUE)

# Generate summary
CST_Summary(data)

## End(Not run)

Function for Calculating the Cluster analysis

Description

This function computes the weather regimes from a cluster analysis. It is applied on the array data in a 's2dv_cube' object. The dimensionality of this object can be also reduced by using PCs obtained from the application of the #'EOFs analysis to filter the dataset. The cluster analysis can be performed with the traditional k-means or those methods included in the hclust (stats package).

Usage

CST_WeatherRegimes(
  data,
  ncenters = NULL,
  EOFs = TRUE,
  neofs = 30,
  varThreshold = NULL,
  method = "kmeans",
  iter.max = 100,
  nstart = 30,
  ncores = NULL
)

Arguments

Value

A list with two elements $data (a 's2dv_cube' object containing the composites cluster = 1,..,K for case (*1) or only k = 1 for any specific cluster, i.e., case (*2)) and $statistics that includes $pvalue (array with the same structure as $data containing the pvalue of the composites obtained through a t-test that accounts for the serial dependence.), cluster (A matrix or vector with integers (from 1:k) indicating the cluster to which each time step is allocated.), persistence (Percentage of days in a month/season before a cluster is replaced for a new one (only if method=’kmeans’ has been selected.)), frequency (Percentage of days in a month/season belonging to each cluster (only if method=’kmeans’ has been selected).),

Author(s)

Verónica Torralba - BSC, veronica.torralba@bsc.es

References

Cortesi, N., V., Torralba, N., González-Reviriego, A., Soret, and F.J., Doblas-Reyes (2019). Characterization of European wind speed variability using weather regimes. Climate Dynamics,53, 4961–4976, doi:10.1007/s00382-019-04839-5.

Torralba, V. (2019) Seasonal climate prediction for the wind energy sector: methods and tools for the development of a climate service. Thesis. Available online: https://eprints.ucm.es/56841/.

Examples

data <- array(abs(rnorm(1280, 283.7, 6)), dim = c(dataset = 2, member = 2, 
                                                 sdate = 3, ftime = 3, 
                                                 lat = 4, lon = 4))
coords <- list(lon = seq(0, 3), lat = seq(47, 44))
obs <- list(data = data, coords = coords)
class(obs) <- 's2dv_cube'

res1 <- CST_WeatherRegimes(data = obs, EOFs = FALSE, ncenters = 4)
res2 <- CST_WeatherRegimes(data = obs, EOFs = TRUE, ncenters = 3)

Forecast Calibration

Description

Five types of member-by-member bias correction can be performed. The "bias" method corrects the bias only, the "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 "mse_min" and "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 "mse_min" method minimizes a constrained mean-squared error using three parameters, the "crps_min" method features four parameters and minimizes the Continuous Ranked Probability Score (CRPS). The "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-k-out cross validation) calibration are possible.

Usage

Calibration(
  exp,
  obs,
  exp_cor = NULL,
  cal.method = "mse_min",
  eval.method = "leave-k-out",
  multi.model = FALSE,
  na.fill = TRUE,
  na.rm = TRUE,
  apply_to = NULL,
  alpha = NULL,
  memb_dim = "member",
  sdate_dim = "sdate",
  dat_dim = NULL,
  ncores = NULL,
  k = 1,
  tail.out = TRUE
)

Arguments

Details

Both the na.fill and na.rm parameters can be used to indicate how the function has to handle the NA values. The na.fill parameter checks whether there are more than three forecast-observations pairs to perform the computation. In case there are three or less pairs, the computation is not carried out, and the value returned by the function depends on the value of this parameter (either NA if na.fill == TRUE or the uncorrected value if na.fill == TRUE). On the other hand, na.rm is used to indicate the function whether to remove the missing values during the computation of the parameters needed to perform the calibration.

Value

An array containing the calibrated forecasts with the dimensions nexp, nobs and same dimensions as in the 'exp' array. nexp is the number of experiment (i.e., 'dat_dim' in exp), and nobs is the number of observation (i.e., 'dat_dim' in obs). If dat_dim is NULL, nexp and nobs are omitted. If 'exp_cor' is provided the returned array will be with the same dimensions as 'exp_cor'.

Author(s)

Verónica Torralba, veronica.torralba@bsc.es

Bert Van Schaeybroeck, bertvs@meteo.be

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

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

Van Schaeybroeck, B., & Vannitsem, S. (2011). Post-processing through linear regression. Nonlinear Processes in Geophysics, 18(2), 147. doi:10.5194/npg-18-147-2011

Van Schaeybroeck, B., & Vannitsem, S. (2015). Ensemble post-processing using member-by-member approaches: theoretical aspects. Quarterly Journal of the Royal Meteorological Society, 141(688), 807-818. doi:10.1002/qj.2397

See Also

CST_Start

Examples

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 <- Calibration(exp = mod1, obs = obs1)

Make categorical forecast based on a multi-model forecast with potential for calibrate

Description

This function converts a multi-model ensemble forecast into a categorical forecast by giving the probability for each category. Different methods are available to combine the different ensemble forecasting models into probabilistic categorical forecasts.

See details in ?CST_CategoricalEnsCombination

Usage

CategoricalEnsCombination(
  fc,
  obs,
  cat.method,
  eval.method,
  amt.cat,
  k,
  tail.out,
  ...
)

Arguments

Value

An array containing the categorical forecasts in the element called $data. The first two dimensions of the returned object are named dataset and member and are both of size one. An additional dimension named category is introduced and is of size amt.cat.

Author(s)

Bert Van Schaeybroeck, bertvs@meteo.be

References

Rajagopalan, B., Lall, U., & Zebiak, S. E. (2002). Categorical climate forecasts through regularization and optimal combination of multiple GCM ensembles. Monthly Weather Review, 130(7), 1792-1811.

Robertson, A. W., Lall, U., Zebiak, S. E., & Goddard, L. (2004). Improved combination of multiple atmospheric GCM ensembles for seasonal prediction. Monthly Weather Review, 132(12), 2732-2744.

Van Schaeybroeck, B., & Vannitsem, S. (2019). Postprocessing of Long-Range Forecasts. In Statistical Postprocessing of Ensemble Forecasts (pp. 267-290).

Performing a Bias Correction conditioned by the dynamical properties of the data.

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'

Usage

DynBiasCorrection(
  exp,
  obs,
  method = "QUANT",
  wetday = FALSE,
  proxy = "dim",
  quanti,
  ncores = NULL
)

Arguments

Value

A multidimensional array with named dimensions with a bias correction performed conditioned by local dimension 'dim' or inverse of persistence 'theta'.

Author(s)

Carmen Alvarez-Castro, carmen.alvarez-castro@cmcc.it

Maria M. Chaves-Montero, mdm.chaves-montero@cmcc.it

Veronica Torralba, veronica.torralba@cmcc.it

Davide Faranda, davide.faranda@lsce.ipsl.fr

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: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.

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)

Ensemble clustering

Description

This function performs a clustering on members/starting dates and returns a number of scenarios, with representative members for each of them. The clustering is performed in a reduced EOF space.

Usage

EnsClustering(
  data,
  lat,
  lon,
  time_moment = "mean",
  numclus = NULL,
  lon_lim = NULL,
  lat_lim = NULL,
  variance_explained = 80,
  numpcs = NULL,
  time_percentile = 90,
  time_dim = NULL,
  cluster_dim = "member",
  verbose = TRUE
)

Arguments

Value

A list with elements $cluster (cluster assigned for each member), $freq (relative frequency of each cluster), $closest_member (representative member for each cluster), $repr_field (list of fields for each representative member), composites (list of mean fields for each cluster), $lon (selected longitudes of output fields), $lat (selected longitudes of output fields).

Author(s)

Federico Fabiano - ISAC-CNR, f.fabiano@isac.cnr.it

Ignazio Giuntoli - ISAC-CNR, i.giuntoli@isac.cnr.it

Danila Volpi - ISAC-CNR, d.volpi@isac.cnr.it

Paolo Davini - ISAC-CNR, p.davini@isac.cnr.it

Jost von Hardenberg - ISAC-CNR, j.vonhardenberg@isac.cnr.it

Examples

exp <- array(abs(rnorm(1152))*275, dim = c(dataset = 1, member = 4, 
                                          sdate = 6, ftime = 3, 
                                          lat = 4, lon = 4))
lon <- seq(0, 3)
lat <- seq(48, 45)
res <- EnsClustering(exp, lat = lat, lon = lon, numclus = 2,
                    cluster_dim = c("member", "dataset", "sdate"))

Generate Training and Evaluation Indices for Cross-Validation

Description

This function generates training and evaluation indices based on different cross-validation methods.

Usage

EvalTrainIndices(
  eval.method,
  sample.length,
  sample.length_cor,
  k = 1,
  tail.out = TRUE
)

Arguments

Value

A list of lists, where each element contains: eval.dexes: Indices of evaluation points. train.dexes: Indices of training points.

Author(s)

Theertha Kariyathan, theertha.kariyathan@bsc.es

Examples

# Leave-k-out cross-validation
EvalTrainIndices("leave-k-out", sample.length = 10, sample.length_cor = 5, k = 3)
# Retrospective cross-validation
EvalTrainIndices("retrospective", sample.length = 10, sample.length_cor = 5, k = 3)
# In-sample validation
EvalTrainIndices("in-sample", sample.length = 10, sample.length_cor = 5)
# Hindcast vs. Forecast validation
EvalTrainIndices("hindcast-vs-forecast", sample.length = 10, sample.length_cor = 5)

Function to Split Dimension

Description

This function merges two dimensions of an array into one. The user can select the dimensions to merge and provide the final name of the dimension. The user can select to remove NA values or keep them.

Usage

MergeDims(
  data,
  merge_dims = c("time", "monthly"),
  rename_dim = NULL,
  na.rm = FALSE
)

Arguments

Value

An array with the two selected dimensions merged into a single dimension.

Author(s)

Nuria Perez-Zanon, nuria.perez@bsc.es

Examples

data <- 1 : 20
dim(data) <- c(time = 10, lat = 2)
new_data <- MergeDims(data, merge_dims = c('time', 'lat'))

EOF analysis of multiple variables starting from an array (reduced version)

Description

This function performs EOF analysis over multiple variables, accepting in input an array with a dimension "var" for each variable to analyse. Based on Singular Value Decomposition. For each field the EOFs are computed and the corresponding PCs are standardized (unit variance, zero mean); the minimum number of principal components needed to reach the user-defined variance is retained. The function weights the input data for the latitude cosine square root.

Usage

MultiEOF(
  data,
  lon,
  lat,
  dates,
  time = NULL,
  lon_dim = "lon",
  lat_dim = "lat",
  time_dim = "ftime",
  sdate_dim = "sdate",
  var_dim = "var",
  neof_max = 40,
  neof_composed = 5,
  minvar = 0.6,
  lon_lim = NULL,
  lat_lim = NULL,
  ncores = NULL
)

Arguments

Value

A list containing:

Author(s)

Jost von Hardenberg - ISAC-CNR, j.vonhardenberg@isac.cnr.it

Paolo Davini - ISAC-CNR, p.davini@isac.cnr.it

Examples

exp <- array(runif(1280)*280, dim = c(dataset = 2, member = 2, sdate = 3, 
                                     ftime = 3, lat = 4, lon = 4, var = 1))
lon <- seq(0, 3)
lat <- seq(47, 44)
dates <- c("2000-11-01", "2000-12-01", "2001-01-01", "2001-11-01", 
          "2001-12-01", "2002-01-01", "2002-11-01", "2002-12-01", "2003-01-01")
Dates <- as.POSIXct(dates, format = "%Y-%m-%d")
dim(Dates) <- c(ftime = 3, sdate = 3)
cal <- MultiEOF(data = exp, lon = lon, lat = lat, dates = Dates)

Multiple Metrics applied in Multiple Model Anomalies

Description

This function calculates correlation (Anomaly Correlation Coefficient; ACC), root mean square error (RMS) and the root mean square error skill score (RMSSS) of individual anomaly models and multi-models mean (if desired) with the observations on arrays with named dimensions.

Usage

MultiMetric(
  exp,
  obs,
  metric = "correlation",
  multimodel = TRUE,
  time_dim = "ftime",
  memb_dim = "member",
  sdate_dim = "sdate"
)

Arguments

Value

A list of arrays containing the statistics of the selected metric in the element $data which is a list of arrays: for the metric requested and others for statistics about its signeificance. The arrays have two dataset dimensions equal to the 'dataset' dimension in the exp$data and obs$data inputs. If multimodel is TRUE, the greatest position in the first dimension correspons to the Multi-Model Mean.

Author(s)

Mishra Niti, niti.mishra@bsc.es

Perez-Zanon Nuria, nuria.perez@bsc.es

References

Mishra, N., Prodhomme, C., & Guemas, V. (n.d.). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29-31. doi:10.1007/s00382-018-4404-z

See Also

Corr, RMS, RMSSS and CST_Load

Examples

exp <- array(rnorm(2*2*4*5*2*2), 
            dim = c(dataset = 2, member = 2, sdate = 4, ftime = 5, lat = 2, 
                    lon = 2))
obs <- array(rnorm(1*1*4*5*2*2),
            dim = c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 2, 
                    lon = 2))
res <- MultiMetric(exp = exp, obs = obs)

Computing the Index PDFs for a dataset of SFSs for a hindcats period.

Description

This function implements the computation to obtain the index PDFs (e.g. NAO index) to improve the index estimate from SFSs for a hindcast period.

Usage

PDFIndexHind(
  index_hind,
  index_obs,
  method = "ME",
  time_dim_name = "time",
  na.rm = FALSE
)

Arguments

Value

An array with at least two dimensions (time, statistic = 4). The firt statistic is the parameter 'mean' of the PDF with not bias corrected. The second statistic is the parameter 'standard deviation' of the PDF with not bias corrected. The third statistic is the parameter 'mean' of the PDF with bias corrected. The fourth statistic is the parameter 'standard deviation' of the PDF with bias corrected.

Author(s)

Eroteida Sanchez-Garcia - AEMET, esanchezg@aemet.es

References

Regionally improved seasonal forecast of precipitation through Best estimation of winter NAO, Sanchez-Garcia, E. et al., Adv. Sci. Res., 16, 165174, 2019, doi:10.5194/asr-16-165-2019

Examples

# Example for the PDFIndexHind function
# Example 1 
index_obs <- 1 : (5 * 3 ) 
dim(index_obs) <- c(time = 5, season = 3)
index_hind <- 1 : (5 * 4 * 3)
dim(index_hind) <- c(time = 5, member = 4, season = 3)
res <- PDFIndexHind(index_hind, index_obs)
dim(res)
# time statistic  season
#   5         4        3
# Example 2
index_obs <- 1 : (5 * 3) 
dim(index_obs) <- c(time = 5, season = 3)
index_hind <- 1 : (5 * 2 * 3)
dim(index_hind) <- c(time = 5, statistic = 2, season = 3)
res <- PDFIndexHind(index_hind, index_obs)

Plot Multiple Lon-Lat Variables In a Single Map According to a Decision Function

Description

Plot a number a two dimensional matrices with (longitude, latitude) dimensions on a single map with the cylindrical equidistant latitude and longitude projection.

Usage

PlotCombinedMap(
  maps,
  lon,
  lat,
  map_select_fun,
  display_range,
  map_dim = "map",
  brks = NULL,
  cols = NULL,
  bar_limits = NULL,
  triangle_ends = c(FALSE, FALSE),
  col_inf = NULL,
  col_sup = NULL,
  col_unknown_map = "white",
  mask = NULL,
  col_mask = "grey",
  dots = NULL,
  bar_titles = NULL,
  legend_scale = 1,
  cex_bar_titles = 1.5,
  plot_margin = NULL,
  bar_extra_margin = c(2, 0, 2, 0),
  fileout = NULL,
  width = 8,
  height = 5,
  size_units = "in",
  res = 100,
  drawleg = TRUE,
  return_leg = FALSE,
  ...
)

Arguments

Value

If return_leg = FALSE (default), the function plots a combined map to the current graphical device (or saves it to file if fileout is specified) and returns NULL.

If return_leg = TRUE, no map is plotted. Instead, a list is returned containing the arguments required to plot the color bars manually using s2dv::GradientCatsColorBar or s2dv::ColorBar. This allows users to customize the color bars independently.

Author(s)

Nicolau Manubens, nicolau.manubens@bsc.es

Veronica Torralba, veronica.torralba@bsc.es

See Also

PlotCombinedMap and PlotEquiMap

Examples

# Simple example
x <- array(1:(20 * 10), dim = c(lat = 10, lon = 20)) / 200
a <- x * 0.6
b <- (1 - x) * 0.6
c <- 1 - (a + b)
lons <- seq(0, 359.5, length = 20)
lats <- seq(-89.5, 89.5, length = 10)
PlotCombinedMap(list(a, b, c), lons, lats, 
               toptitle = 'Maximum map',
               map_select_fun = max,
               display_range = c(0, 1),
               bar_titles = paste('% of belonging to', c('a', 'b', 'c')), 
               brks = 20, width = 12, height = 10, 
               fileout = file.path(tempdir(), "example1.pdf"))
unlink(file.path(tempdir(), "example1.pdf"))

Lon <- c(0:40, 350:359)
Lat <- 51:26
data <- rnorm(51 * 26 * 3)
dim(data) <- c(map = 3, lon = 51, lat = 26)
mask <-  sample(c(0,1), replace = TRUE, size = 51 * 26)
dim(mask) <- c(lat = 26, lon = 51)
PlotCombinedMap(data, lon = Lon, lat = Lat, map_select_fun = max,
               display_range = range(data), mask = mask,
               width = 14, height = 10,
               fileout = file.path(tempdir(), "example2.pdf")) 
unlink(file.path(tempdir(), "example2.pdf")) 

Plot one or multiple ensemble forecast pdfs for the same event

Description

This function plots the probability distribution function of several ensemble forecasts. Separate panels are used to plot forecasts valid or initialized at different times or by different models or even at different locations. Probabilities for tercile categories are computed, plotted in colors and annotated. An asterisk marks the tercile with higher probabilities. Probabilities for extreme categories (above P90 and below P10) can also be included as hatched areas. Individual ensemble members can be plotted as jittered points. The observed value is optionally shown as a diamond.

Usage

PlotForecastPDF(
  fcst,
  tercile.limits,
  extreme.limits = NULL,
  obs = NULL,
  plotfile = NULL,
  title = "Set a title",
  var.name = "Varname (units)",
  fcst.names = NULL,
  add.ensmemb = c("above", "below", "no"),
  color.set = c("ggplot", "s2s4e", "hydro", "vitigeoss"),
  memb_dim = "member"
)

Arguments

Value

A ggplot object containing the plot.

Author(s)

Llorenç Lledó llledo@bsc.es

Examples

fcsts <- data.frame(fcst1 = rnorm(10), fcst2 = rnorm(10, 0.5, 1.2))
PlotForecastPDF(fcsts,c(-1,1))

Plot Maps of Most Likely Quantiles

Description

This function receives as main input (via the parameter probs) a collection of longitude-latitude maps, each containing the probabilities (from 0 to 1) of the different grid cells of belonging to a category. As many categories as maps provided as inputs are understood to exist. The maps of probabilities must be provided on a common rectangular regular grid, and a vector with the longitudes and a vector with the latitudes of the grid must be provided. The input maps can be provided in two forms, either as a list of multiple two-dimensional arrays (one for each category) or as a three-dimensional array, where one of the dimensions corresponds to the different categories.

Usage

PlotMostLikelyQuantileMap(
  probs,
  lon,
  lat,
  cat_dim = "bin",
  bar_titles = NULL,
  col_unknown_cat = "white",
  drawleg = TRUE,
  ...
)

Arguments

Value

Produces a map showing the most likely category or quantile for each grid cell, based on the input probability maps.

Author(s)

Veronica Torralba, veronica.torralba@bsc.es, Nicolau Manubens, nicolau.manubens@bsc.es

See Also

PlotCombinedMap and PlotEquiMap

Examples

# Simple example
x <- array(1:(20 * 10), dim = c(lat = 10, lon = 20)) / 200
a <- x * 0.6
b <- (1 - x) * 0.6
c <- 1 - (a + b)
lons <- seq(0, 359.5, length = 20)
lats <- seq(-89.5, 89.5, length = 10)
PlotMostLikelyQuantileMap(list(a, b, c), lons, lats, 
                         toptitle = 'Most likely tercile map',
                         bar_titles = paste('% of belonging to', c('a', 'b', 'c')), 
                         brks = 20, width = 10, height = 8,
                         bar_extra_margin = c(1, 1, 1, 1),
                         fileout = file.path(tempdir(),"example1.pdf"))
unlink(file.path(tempdir(),"example1.pdf"))

# More complex example
n_lons <- 40
n_lats <- 20
n_timesteps <- 100
n_bins <- 4

# 1. Generation of sample data
lons <- seq(0, 359.5, length = n_lons)
lats <- seq(-89.5, 89.5, length = n_lats)

# This function builds a 3-D gaussian at a specified point in the map.
make_gaussian <- function(lon, sd_lon, lat, sd_lat) {
 w <- outer(lons, lats, function(x, y) dnorm(x, lon, sd_lon) * dnorm(y, lat, sd_lat))
 min_w <- min(w)
 w <- w - min_w
 w <- w / max(w)
 w <- t(w)
 names(dim(w)) <- c('lat', 'lon')
 w
}

# This function generates random time series (with values ranging 1 to 5) 
# according to 2 input weights.
gen_data <- function(w1, w2, n) {
 r <- sample(1:5, n, 
             prob = c(.05, .9 * w1, .05, .05, .9 * w2), 
             replace = TRUE)
 r <- r + runif(n, -0.5, 0.5)
 dim(r) <- c(time = n)
 r
}

# We build two 3-D gaussians.
w1 <- make_gaussian(120, 80, 20, 30)
w2 <- make_gaussian(260, 60, -10, 40)

# We generate sample data (with dimensions time, lat, lon) according 
# to the generated gaussians
sample_data <- multiApply::Apply(list(w1, w2), NULL, 
                                gen_data, n = n_timesteps)$output1

# 2. Binning sample data
prob_thresholds <- 1:n_bins / n_bins
prob_thresholds <- prob_thresholds[1:(n_bins - 1)]
thresholds <- quantile(sample_data, prob_thresholds)

binning <- function(x, thresholds) {
 n_samples <- length(x)
 n_bins <- length(thresholds) + 1

 thresholds <- c(thresholds, max(x))
 result <- 1:n_bins
 lower_threshold <- min(x) - 1
 for (i in 1:n_bins) {
   result[i] <- sum(x > lower_threshold & x <= thresholds[i]) / n_samples
   lower_threshold <- thresholds[i]
 }

 dim(result) <- c(bin = n_bins)
 result
}

bins <- multiApply::Apply(sample_data, 'time', binning, thresholds)$output1

# 3. Plotting most likely quantile/bin
PlotMostLikelyQuantileMap(bins, lons, lats, 
                         toptitle = 'Most likely quantile map',
                         bar_titles = paste('% of belonging to', letters[1:n_bins]),
                         mask = 1 - (w1 + w2 / max(c(w1, w2))), 
                         brks = 20, width = 10, height = 8,
                         bar_extra_margin = c(1, 1, 1, 1),
                         fileout = file.path(tempdir(),"example2.pdf"))
unlink(file.path(tempdir(),"example2.pdf"))

Plotting two probability density gaussian functions and the optimal linear estimation (OLE) as result of combining them.

Description

This function plots two probability density gaussian functions and the optimal linear estimation (OLE) as result of combining them.

Usage

PlotPDFsOLE(
  pdf_1,
  pdf_2,
  nsigma = 3,
  legendPos = "bottom",
  legendSize = 1,
  plotfile = NULL,
  width = 30,
  height = 15,
  units = "cm",
  dpi = 300
)

Arguments

Value

PlotPDFsOLE() returns a ggplot object containing the plot.

Author(s)

Eroteida Sanchez-Garcia - AEMET, esanchezg@aemet.es

Examples

# Example 1
pdf_1 <- c(1.1,0.6)
attr(pdf_1, "name") <- "NAO1"
dim(pdf_1) <-  c(statistic = 2)
pdf_2 <- c(1,0.5)
attr(pdf_2, "name") <- "NAO2"
dim(pdf_2) <-  c(statistic = 2)

PlotPDFsOLE(pdf_1, pdf_2)

Function to convert any 3-d numerical array to a grid of coloured triangles.

Description

This function converts a 3-d numerical data array into a coloured grid with triangles. It is useful for a slide or article to present tabular results as colors instead of numbers. This can be used to compare the outputs of two or four categories (e.g. modes of variability, clusters, or forecast systems).

Usage

PlotTriangles4Categories(
  data,
  brks = NULL,
  cols = NULL,
  toptitle = NULL,
  sig_data = NULL,
  pch_sig = 18,
  col_sig = "black",
  cex_sig = 1,
  xlab = TRUE,
  ylab = TRUE,
  xlabels = NULL,
  xtitle = NULL,
  ylabels = NULL,
  ytitle = NULL,
  legend = TRUE,
  lab_legend = NULL,
  cex_leg = 1,
  col_leg = "black",
  cex_axis = 1.5,
  mar = c(5, 4, 0, 0),
  fileout = NULL,
  size_units = "px",
  res = 100,
  figure.width = 1,
  ...
)

Arguments

Value

A figure in popup window by default, or saved to the specified path.

Author(s)

History:
1.0 - 2020-10 (V.Torralba, veronica.torralba@bsc.es) - Original code

Examples

# Example with random data
arr1 <- array(runif(n = 4 * 5 * 4, min = -1, max = 1), dim = c(4,5,4))
names(dim(arr1)) <- c('dimx', 'dimy', 'dimcat')
arr2 <- array(TRUE, dim = dim(arr1))
arr2[which(arr1 < 0.3)] <- FALSE
PlotTriangles4Categories(data = arr1,
                        cols = c('white','#fef0d9','#fdd49e','#fdbb84','#fc8d59'),
                        brks = c(-1, 0, 0.1, 0.2, 0.3, 0.4), 
                        lab_legend = c('NAO+', 'BL','AR','NAO-'), 
                        xtitle = "Target month", ytitle = "Lead time",
                        xlabels = c("Jan", "Feb", "Mar", "Apr"))

Plots the observed weekly means and climatology of a timeseries data

Description

This function plots the observed weekly means and climatology of a timeseries data using ggplot package. It compares the weekly climatology in a specified period (reference period) to the observed conditions during the target period analyzed in the case study.

Usage

PlotWeeklyClim(
  data,
  first_date,
  ref_period,
  last_date = NULL,
  data_years = NULL,
  time_dim = "time",
  sdate_dim = "sdate",
  ylim = NULL,
  title = NULL,
  subtitle = NULL,
  ytitle = NULL,
  legend = TRUE,
  palette = "Blues",
  fileout = NULL,
  device = NULL,
  width = 8,
  height = 6,
  units = "in",
  dpi = 300
)

Arguments

Value

A ggplot object containing the plot.

Examples

data <- array(rnorm(49*20*3, 274), dim = c(time = 49, sdate = 20, member = 3))
PlotWeeklyClim(data = data, first_date = '2002-08-09', 
              last_date = '2002-09-15', ref_period = 2010:2019, 
              data_years = 2000:2019, time_dim = 'time', sdate_dim = 'sdate',
              title = "Observed weekly means and climatology", 
              subtitle = "Target years: 2010 to 2019", 
              ytitle = paste0('tas', " (", "deg.C", ")"))

Computing scores of predictability using two dynamical proxies based on dynamical systems theory.

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.

Usage

Predictability(dim, theta, ncores = NULL)

Arguments

Value

A list of length 2:

• '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'.

• '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.

Author(s)

Carmen Alvarez-Castro, carmen.alvarez-castro@cmcc.it

Maria M. Chaves-Montero, mdm.chaves-montero@cmcc.it

Veronica Torralba, veronica.torralba@cmcc.it

Davide Faranda, davide.faranda@lsce.ipsl.fr

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: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.

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)

Computing two dinamical proxies of the attractor.

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:

Usage

ProxiesAttractor(data, quanti, ncores = NULL)

Arguments

Value

dim and theta

Author(s)

Carmen Alvarez-Castro, carmen.alvarez-castro@cmcc.it

Maria M. Chaves-Montero, mdm.chaves-montero@cmcc.it

Veronica Torralba, veronica.torralba@cmcc.it

Davide Faranda, davide.faranda@lsce.ipsl.fr

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: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.

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))
plot(time, Attractor$dim, xlab = 'time', ylab = 'd',
    main = 'local dimension', type = 'l')

Quantile Mapping for seasonal or decadal forecast data

Description

This function is a wrapper of fitQmap and doQmap from package 'qmap' to be applied on multi-dimensional arrays. The quantile mapping adjustment between an experiment, typically a hindcast, and observation is applied to the experiment itself or to a provided forecast.

Usage

QuantileMapping(
  exp,
  obs,
  exp_cor = NULL,
  sdate_dim = "sdate",
  memb_dim = "member",
  window_dim = NULL,
  method = "QUANT",
  na.rm = FALSE,
  eval.method = "leave-k-out",
  k = 1,
  ncores = NULL,
  ...
)

Arguments

Value

An array containing the experimental data after applying the quantile mapping correction.

Author(s)

Nuria Perez-Zanon, nuria.perez@bsc.es

See Also

fitQmap and doQmap

Examples

# Use synthetic data
set.seed(123)
exp <- as.numeric(1:prod(6,10,15))
dim(exp) <- c(member = 6, syear = 10, window = 15)
obs <- as.numeric(rnorm(prod(1,10,15), 50))
dim(obs) <- c(member = 1, syear = 10, window = 15)
fcst <- 100*(1:prod(8,1,1))
dim(fcst) <- c(member = 8, syear = 1, swindow = 1)
res <- QuantileMapping(exp = exp, obs = obs, exp_cor = fcst, 
                      memb_dim = 'member', sdate_dim = 'syear', window_dim = 'window')

RainFARM spectral slopes from an array (reduced version)

Description

This function computes spatial spectral slopes from an array, to be used for RainFARM stochastic precipitation downscaling method.

Usage

RFSlope(
  data,
  kmin = 1,
  time_dim = NULL,
  lon_dim = "lon",
  lat_dim = "lat",
  ncores = NULL
)

Arguments

Value

RFSlope() returns spectral slopes using the RainFARM convention (the logarithmic slope of k*|A(k)|^2 where A(k) are the spectral amplitudes). The returned array has the same dimensions as the input array, minus the dimensions specified by lon_dim, lat_dim and time_dim.

Author(s)

Jost von Hardenberg - ISAC-CNR, j.vonhardenberg@isac.cnr.it

Examples

# Example for the 'reduced' RFSlope function 
# Create a test array with dimension 8x8 and 20 timesteps, 
# 3 starting dates and 20 ensemble members.
pr <- 1:(4*3*8*8*20)
dim(pr) <- c(ensemble = 4, sdate = 3, lon = 8, lat = 8, ftime = 20)
# Compute the spectral slopes ignoring the wavenumber
# corresponding to the largest scale (the box)
slopes <- RFSlope(pr, kmin = 2, time_dim = 'ftime')

Temperature downscaling of a CSTools object using lapse rate correction (reduced version)

Description

This function implements a simple lapse rate correction of a temperature field (a multidimensional array) as input. The input lon grid must be increasing (but can be modulo 360). The input lat grid can be irregularly spaced (e.g. a Gaussian grid) The output grid can be irregularly spaced in lon and/or lat.

Usage

RFTemp(
  data,
  lon,
  lat,
  oro,
  lonoro,
  latoro,
  xlim = NULL,
  ylim = NULL,
  lapse = 6.5,
  lon_dim = "lon",
  lat_dim = "lat",
  time_dim = NULL,
  nolapse = FALSE,
  verbose = FALSE,
  compute_delta = FALSE,
  method = "bilinear",
  delta = NULL
)

Arguments

Value

CST_RFTemp() returns a downscaled CSTools object.

RFTemp() returns a list containing the fine-scale longitudes, latitudes and the downscaled fields.

Author(s)

Jost von Hardenberg - ISAC-CNR, j.vonhardenberg@isac.cnr.it

References

Method described in ERA4CS MEDSCOPE milestone M3.2: High-quality climate prediction data available to WP4 here: https://www.medscope-project.eu/the-project/deliverables-reports/ and in H2020 ECOPOTENTIAL Deliverable No. 8.1: High resolution (1-10 km) climate, land use and ocean change scenarios here: https://ec.europa.eu/research/participants/documents/downloadPublic?documentIds=080166e5b6cd2324&appId=PPGMS.

Examples

# Generate simple synthetic data and downscale by factor 4
t <- rnorm(7 * 6 * 4 * 3) * 10 + 273.15 + 10
dim(t) <- c(sdate = 3, ftime = 4, lat = 6, lon = 7)
lon <- seq(3, 9, 1)
lat <- seq(42, 47, 1)
o <- runif(29 * 29) * 3000
dim(o) <- c(lat = 29, lon = 29)
lono <- seq(3, 10, 0.25)
lato <- seq(41, 48, 0.25)
res <- RFTemp(t, lon, lat, o, lono, lato, xlim = c(4, 8), ylim = c(43, 46),
             lapse = 6.5, time_dim = 'ftime')

Compute climatological weights for RainFARM stochastic precipitation downscaling

Description

Compute climatological ("orographic") weights from a fine-scale precipitation climatology file.

Usage

RF_Weights(
  zclim,
  latin,
  lonin,
  nf,
  lat,
  lon,
  fsmooth = TRUE,
  lonname = "lon",
  latname = "lat",
  ncores = NULL
)

Arguments

Value

An object of class 's2dv_cube' containing in matrix data the weights with dimensions (lon, lat).

Author(s)

Jost von Hardenberg - ISAC-CNR, j.vonhardenberg@isac.cnr.it

References

Terzago, S., Palazzi, E., & von Hardenberg, J. (2018). Stochastic downscaling of precipitation in complex orography: A simple method to reproduce a realistic fine-scale climatology. Natural Hazards and Earth System Sciences, 18(11), 2825-2840. doi:10.5194/nhess-18-2825-2018.

Examples

a <- array(1:2500, c(lat = 50, lon = 50))
res <- RF_Weights(a, seq(0.1 ,5, 0.1), seq(0.1 ,5, 0.1), 
                 nf = 5, lat = 1:5, lon = 1:5)

RainFARM stochastic precipitation downscaling (reduced version)

Description

This function implements the RainFARM stochastic precipitation downscaling method and accepts in input an array with named dims ("lon", "lat") and one or more dimension (such as "ftime", "sdate" or "time") over which to average automatically determined spectral slopes. Adapted for climate downscaling and including orographic correction. References: Terzago, S. et al. (2018). NHESS 18(11), 2825-2840. doi:10.5194/nhess-18-2825-2018, D'Onofrio et al. (2014), J of Hydrometeorology 15, 830-843; Rebora et. al. (2006), JHM 7, 724.

Usage

RainFARM(
  data,
  lon,
  lat,
  nf,
  weights = 1,
  nens = 1,
  slope = 0,
  kmin = 1,
  fglob = FALSE,
  fsmooth = TRUE,
  nprocs = 1,
  time_dim = NULL,
  lon_dim = "lon",
  lat_dim = "lat",
  drop_realization_dim = FALSE,
  verbose = FALSE
)

Arguments

Details

Wether parameter 'slope' and 'weights' presents seasonality dependency, a dimension name should match between these parameters and the input data in parameter 'data'. See example 2 below where weights and slope vary with 'sdate' dimension.

Value

RainFARM() Returns a list containing the fine-scale longitudes, latitudes and the sequence of nens downscaled fields. If nens > 1 an additional dimension named "realization" is added to the output array after the "member" dimension (if it exists and unless drop_realization_dim = TRUE is specified). The ordering of the remaining dimensions in the exp element of the input object is maintained.

Author(s)

Jost von Hardenberg - ISAC-CNR, j.vonhardenberg@isac.cnr.it

Examples

# Example for the 'reduced' RainFARM function 
nf <- 8   # Choose a downscaling by factor 8
exp <- 1 : (2 * 3 * 4 * 8 * 8)
dim(exp) <- c(dataset = 1, member = 2, sdate = 3, ftime = 4, lat = 8, lon = 8)
lon <- seq(10, 13.5, 0.5)
lat <- seq(40, 43.5, 0.5)
# Create a test array of weights
ww <- array(1., dim = c(lon = 8 * nf, lat = 8 * nf))
res <- RainFARM(data = exp, lon = lon, lat = lat, nf = nf, 
               weights = ww, nens = 3, time_dim = 'ftime')

Function for matching a field of anomalies with a set of maps used as a reference (e.g. clusters obtained from the WeatherRegime function).

Description

This function performs the matching between a field of anomalies and a set of maps which will be used as a reference. The anomalies will be assigned to the reference map for which the minimum Eucledian distance (method = 'distance') or highest spatial correlation (method = 'ACC') is obtained.

Usage

RegimesAssign(
  data,
  ref_maps,
  lat,
  method = "distance",
  composite = FALSE,
  memb = FALSE,
  ncores = NULL
)

Arguments

Value

A list with elements $composite (3-d array (lon, lat, k) containing the composites k = 1,..,K for case (*1) or only k = 1 for any specific cluster, i.e., case (*2)) (only if composite = 'TRUE'), $pvalue (array with the same structure as $composite containing the pvalue of the composites obtained through a t-test that accounts for the serial dependence of the data with the same structure as Composite.) (only if composite='TRUE'), $cluster (array with the same dimensions as data (except latitude and longitude which are removed) indicating the ref_maps to which each point is allocated.), $frequency (A vector of integers (from k = 1, ... k n reference maps) indicating the percentage of assignations corresponding to each map.),

Author(s)

Verónica Torralba - BSC, veronica.torralba@bsc.es

References

Torralba, V. (2019) Seasonal climate prediction for the wind energy sector: methods and tools for the development of a climate service. Thesis. Available online: https://eprints.ucm.es/56841/

Examples

data <- array(abs(rnorm(1280, 282.7, 6.4)), dim = c(dataset = 2, member = 2, 
                                                   sdate = 3, ftime = 3, 
                                                   lat = 4, lon = 4))
regimes <- WeatherRegime(data = data, lat = seq(47, 44),
                        EOFs = FALSE, ncenters = 4)$composite
res1 <- RegimesAssign(data = data, ref_maps = drop(regimes), 
                     lat = seq(47, 44), composite = FALSE)

Save a multidimensional array with metadata to data in NetCDF format

Description

This function allows to save a data array with metadata into a NetCDF file, allowing to reload the saved data using Start function from StartR package. If the original 's2dv_cube' object has been created from CST_Load(), then it can be reloaded with Load().

Usage

SaveExp(
  data,
  destination = tempdir(),
  coords = NULL,
  Dates = NULL,
  time_bounds = NULL,
  startdates = NULL,
  varname = NULL,
  metadata = NULL,
  Datasets = NULL,
  sdate_dim = "sdate",
  ftime_dim = "time",
  memb_dim = "member",
  dat_dim = "dataset",
  var_dim = "var",
  drop_dims = NULL,
  single_file = FALSE,
  extra_string = NULL,
  global_attrs = NULL,
  units_hours_since = FALSE
)

Arguments

Value

Multiple or single NetCDF files containing the data array.

Author(s)

Perez-Zanon Nuria, nuria.perez@bsc.es

Examples

data <- lonlat_temp_st$exp$data
lon <- lonlat_temp_st$exp$coords$lon
lat <- lonlat_temp_st$exp$coords$lat
coords <- list(lon = lon, lat = lat)
Datasets <- lonlat_temp_st$exp$attrs$Datasets
varname <- 'tas'
Dates <- lonlat_temp_st$exp$attrs$Dates
metadata <- lonlat_temp_st$exp$attrs$Variable$metadata
SaveExp(data = data, coords = coords, Datasets = Datasets, varname = varname, 
       Dates = Dates, metadata = metadata, single_file = TRUE, 
       ftime_dim = 'ftime', var_dim = 'var', dat_dim = 'dataset')
       

Function to Split Dimension

Description

This function split a dimension in two. The user can select the dimension to split and provide indices indicating how to split that dimension or dates and the frequency expected (monthly or by day, month and year). The user can also provide a numeric frequency indicating the length of each division.

Usage

SplitDim(
  data,
  split_dim = "time",
  indices,
  freq = "monthly",
  new_dim_name = NULL,
  dates = NULL,
  return_indices = FALSE
)

Arguments

Value

An n-dimensional array with the specified dimension split into a new dimension, preserving all other original dimensions and their names,

Author(s)

Nuria Perez-Zanon, nuria.perez@bsc.es

Examples

data <- 1 : 20
dim(data) <- c(time = 10, lat = 2)
indices <- c(rep(1,5), rep(2,5))
new_data <- SplitDim(data, indices = indices)
time <- c(seq(ISOdate(1903, 1, 1), ISOdate(1903, 1, 4), "days"),
         seq(ISOdate(1903, 2, 1), ISOdate(1903, 2, 4), "days"),
         seq(ISOdate(1904, 1, 1), ISOdate(1904, 1, 2), "days"))
new_data <- SplitDim(data, indices = time)
new_data <- SplitDim(data, indices = time, freq = 'day')
new_data <- SplitDim(data, indices = time, freq = 'month')
new_data <- SplitDim(data, indices = time, freq = 'year')

Function for Calculating the Cluster analysis

Description

This function computes the weather regimes from a cluster analysis. It can be applied over the dataset with dimensions c(year/month, month/day, lon, lat), or by using PCs obtained from the application of the EOFs analysis to filter the dataset. The cluster analysis can be performed with the traditional k-means or those methods included in the hclust (stats package).

Usage

WeatherRegime(
  data,
  ncenters = NULL,
  EOFs = TRUE,
  neofs = 30,
  varThreshold = NULL,
  lon = NULL,
  lat = NULL,
  method = "kmeans",
  iter.max = 100,
  nstart = 30,
  ncores = NULL
)

Arguments

Value

A list with elements $composite (array with at least 3-d ('lat', 'lon', 'cluster') containing the composites k = 1,..,K for case (*1) or only k = 1 for any specific cluster, i.e., case (*2)), pvalue (array with at least 3-d ('lat','lon','cluster') with the pvalue of the composites obtained through a t-test that accounts for the serial dependence of the data with the same structure as Composite.), cluster (A matrix or vector with integers (from 1:k) indicating the cluster to which each time step is allocated.), persistence (Percentage of days in a month/season before a cluster is replaced for a new one (only if method=’kmeans’ has been selected.)), frequency (Percentage of days in a month/season belonging to each cluster (only if method=’kmeans’ has been selected).),

Author(s)

Verónica Torralba - BSC, veronica.torralba@bsc.es

References

Cortesi, N., V., Torralba, N., González-Reviriego, A., Soret, and F.J., Doblas-Reyes (2019). Characterization of European wind speed variability using weather regimes. Climate Dynamics,53, 4961–4976, doi:10.1007/s00382-019-04839-5.

Torralba, V. (2019) Seasonal climate prediction for the wind energy sector: methods and tools for the development of a climate service. Thesis. Available online: https://eprints.ucm.es/56841/

Examples

data <- array(abs(rnorm(1280, 283.7, 6)), dim = c(dataset = 2, member = 2,  
                                                 sdate = 3, ftime = 3, 
                                                 lat = 4, lon = 4))
lat <- seq(47, 44)
res <- WeatherRegime(data = data, lat = lat,  
                    EOFs = FALSE, ncenters = 4)

Conversion of 'startR_array' or 'list' objects to 's2dv_cube'

Description

This function converts data loaded using Start function from startR package or Load from s2dv into an 's2dv_cube' object.

Usage

as.s2dv_cube(object, remove_attrs_coords = FALSE, remove_null = FALSE)

Arguments

Value

The function returns an 's2dv_cube' object to be easily used with functions with the prefix CST from CSTools and CSIndicators packages. The object is mainly a list with the following elements:

• 'data', array with named dimensions;

• 'dims', named vector of the data dimensions;

• 'coords', list of named vectors with the coordinates corresponding to the dimensions of the data parameter;

• 'attrs', named list with elements:

  • 'Dates', array with named temporal dimensions of class 'POSIXct' from time values in the data;

  • 'Variable', has the following components:

    • 'varName', character vector of the short variable name. It is usually specified in the parameter 'var' from the functions Start and Load;

    • 'metadata', named list of elements with variable metadata. They can be from coordinates variables (e.g. longitude) or main variables (e.g. 'var');

  • 'Datasets', character strings indicating the names of the datasets;

  • 'source_files', a vector of character strings with complete paths to all the found files involved in loading the data;

  • 'when', a time stamp of the date issued by the Start() or Load() call to obtain the data;

  • 'load_parameters', it contains the components used in the arguments to load the data from Start() or Load() functions.

Author(s)

Perez-Zanon Nuria, nuria.perez@bsc.es

Nicolau Manubens, nicolau.manubens@bsc.es

See Also

s2dv_cube, CST_Start, Start and CST_Load

Examples

## Not run: 
# Example 1: convert an object from startR::Start function to 's2dv_cube'
library(startR)
repos <- '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc'
data <- Start(dat = repos,
             var = 'tas',
             sdate = c('20170101', '20180101'),
             ensemble = indices(1:5),
             time = 'all',
             latitude = indices(1:5),
             longitude = indices(1:5),
             return_vars = list(latitude = 'dat', longitude = 'dat', time = 'sdate'),
             retrieve = TRUE)
data <- as.s2dv_cube(data)
# Example 2: convert an object from s2dv::Load function to 's2dv_cube'
startDates <- c('20001101', '20011101', '20021101',
               '20031101', '20041101', '20051101')
data <- Load(var = 'tas', exp = 'system5c3s', 
            nmember = 2, sdates = startDates,
            leadtimemax = 3, latmin = 10, latmax = 30,
            lonmin = -10, lonmax = 10, output = 'lonlat')
data <- as.s2dv_cube(data)

## End(Not run)

Sample Of Experimental Precipitation Data In Function Of Longitudes And Latitudes

Description

This sample data set contains a small cutout of gridded seasonal precipitation forecast data from the Copernicus Climate Change ECMWF-System 5 forecast system, to be used to demonstrate downscaling. Specifically, for the 'pr' (precipitation) variable, for the first 6 forecast ensemble members, daily values, for all 31 days in March following the forecast starting dates in November of years 2010 to 2012, for a small 4x4 pixel cutout in a region in the North-Western Italian Alps (44N-47N, 6E-9E). The data resolution is 1 degree.

Details

The 'CST_Load' call used to generate the data set in the infrastructure of the Marconi machine at CINECA is shown next, working on files which were extracted from forecast data available in the MEDSCOPE internal archive.

Value

sample s2dv_cube object

library(CSTools)
infile <- list(path = paste0('/esarchive/exp/ecmwf/system5c3s/daily_mean/',
                             '$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc'))
lonlat_prec <- CST_Load('prlr', exp = list(infile), obs = NULL,
                        sdates = c('20101101', '20111101', '20121101'),
                        leadtimemin = 121, leadtimemax = 151,
                        latmin = 44, latmax = 47,
                        lonmin = 6, lonmax = 9,
                        nmember = 6,
                        storefreq = "daily", sampleperiod = 1,
                        output = "lonlat"
                       )

Author(s)

Jost von Hardenberg j.vonhardenberg@isac.cnr.it

An-Chi Ho an.ho@bsc.es

Sample Of Experimental Precipitation Data In Function Of Longitudes And Latitudes with Start

Description

This sample data set contains a small cutout of gridded seasonal precipitation forecast data from the Copernicus Climate Change ECMWF-System 5 forecast system, to be used to demonstrate downscaling. Specifically, for the 'pr' (precipitation) variable, for the first 6 forecast ensemble members, daily values, for all 31 days in March following the forecast starting dates in November of years 2010 to 2012, for a small 4x4 pixel cutout in a region in the North-Western Italian Alps (44N-47N, 6E-9E). The data resolution is 1 degree.

Details

The 'CST_Start' call used to generate the data set in the infrastructure of the Marconi machine at CINECA is shown next, working on files which were extracted from forecast data available in the MEDSCOPE internal archive.

Value

sample s2dv_cube object

 path <- paste0('/esarchive/exp/ecmwf/system5c3s/daily_mean/',
                '$var$_s0-24h/$var$_$sdate$.nc')
 sdates = c('20101101', '20111101', '20121101')
 latmin <- 44
 latmax <- 47
 lonmin <- 6
 lonmax <- 9

 lonlat_prec_st <- CST_Start(dataset = path,
                             var = 'prlr',
                             member = startR::indices(1:6),
                             sdate = sdates,
                             ftime = 121:151,
                             lat = startR::values(list(latmin, latmax)),
                             lat_reorder = startR::Sort(decreasing = TRUE),
                             lon = startR::values(list(lonmin, lonmax)),
                             lon_reorder = startR::CircularSort(0, 360),
                             synonims = list(lon = c('lon', 'longitude'),
                                             lat = c('lat', 'latitude'),
                                             ftime = c('time', 'ftime'),
                                             member = c('member', 'ensemble')),
                              return_vars = list(ftime = 'sdate',
                                                 lon = NULL, lat = NULL),
                              retrieve = TRUE)

Author(s)

Jost von Hardenberg j.vonhardenberg@isac.cnr.it

An-Chi Ho an.ho@bsc.es

Sample Of Experimental And Observational Climate Data In Function Of Longitudes And Latitudes

Description

This sample data set contains gridded seasonal forecast and corresponding observational data from the Copernicus Climate Change ECMWF-System 5 forecast system, and from the Copernicus Climate Change ERA-5 reconstruction. Specifically, for the 'tas' (2-meter temperature) variable, for the 15 first forecast ensemble members, monthly averaged, for the 3 first forecast time steps (lead months 1 to 4) of the November start dates of 2000 to 2005, for the Mediterranean region (27N-48N, 12W-40E). The data was generated on (or interpolated onto, for the reconstruction) a rectangular regular grid of size 360 by 181.

Details

It is recommended to use the data set as follows:

Value

sample s2dv_cube object

require(zeallot)
c(exp, obs) 

The 'CST_Load' call used to generate the data set in the infrastructure of the Earth Sciences Department of the Barcelona Supercomputing Center is shown next. Note that 'CST_Load' internally calls 's2dv::Load', which would require a configuration file (not provided here) expressing the distribution of the 'system5c3s' and 'era5' NetCDF files in the file system.

library(CSTools)
require(zeallot)

startDates <- c('20001101', '20011101', '20021101',
                '20031101', '20041101', '20051101')

lonlat_temp <- 
  CST_Load(
    var = 'tas', 
    exp = 'system5c3s', 
    obs = 'era5', 
    nmember = 15,
    sdates = startDates,
    leadtimemax = 3,
    latmin = 27, latmax = 48,
    lonmin = -12, lonmax = 40, 
    output = 'lonlat',
    nprocs = 1
  )

Author(s)

Nicolau Manubens nicolau.manubens@bsc.es

Sample Of Experimental And Observational Climate Data In Function Of Longitudes And Latitudes with Start

Description

This sample data set contains gridded seasonal forecast and corresponding observational data from the Copernicus Climate Change ECMWF-System 5 forecast system, and from the Copernicus Climate Change ERA-5 reconstruction. Specifically, for the 'tas' (2-meter temperature) variable, for the 15 first forecast ensemble members, monthly averaged, for the 3 first forecast time steps (lead months 1 to 4) of the November start dates of 2000 to 2005, for the Mediterranean region (27N-48N, 12W-40E). The data was generated on (or interpolated onto, for the reconstruction) a rectangular regular grid of size 360 by 181.

Details

The 'CST_Start' call used to generate the data set in the infrastructure of the Earth Sciences Department of the Barcelona Supercomputing Center is shown next. Note that 'CST_Start' internally calls 'startR::Start' and then uses 'as.s2dv_cube' that converts the 'startR_array' into 's2dv_cube'.

Value

sample s2dv_cube object

 lonlat_temp_st <- NULL
 repos_exp <- paste0('/esarchive/exp/ecmwf/system5c3s/monthly_mean/',
                     '$var$_f6h/$var$_$sdate$.nc')
 sdates <- sapply(2000:2005, function(x) paste0(x, '1101'))
 lonmax <- 40
 lonmin <- -12
 latmax <- 48
 latmin <- 27
 lonlat_temp_st$exp <- CST_Start(dataset = repos_exp,
                                 var = 'tas',
                                 member = startR::indices(1:15),
                                 sdate = sdates,
                                 ftime = startR::indices(1:3),
                                 lat = startR::values(list(latmin, latmax)),
                                 lat_reorder = startR::Sort(decreasing = TRUE), 
                                 lon = startR::values(list(lonmin, lonmax)),
                                 lon_reorder = startR::CircularSort(0, 360),
                                 synonims = list(lon = c('lon', 'longitude'),
                                                 lat = c('lat', 'latitude'),
                                                 member = c('member', 'ensemble'),
                                                 ftime = c('ftime', 'time')),
                                 return_vars = list(lat = NULL, 
                                                    lon = NULL, ftime = 'sdate'),
                                                    retrieve = TRUE)
 
 dates <- c(paste0(2000, c(11, 12)), paste0(2001, c('01', 11, 12)),
            paste0(2002, c('01', 11, 12)),  paste0(2003, c('01', 11, 12)),
            paste0(2004, c('01', 11, 12)),  paste0(2005, c('01', 11, 12)), 200601)
 dates <- sapply(dates, function(x) {paste0(x, '01')})
 dates <- as.POSIXct(dates, format = '
 dim(dates) <- c(ftime = 3, sdate = 6)

 dates <- t(dates)
 names(dim(dates)) <- c('sdate', 'ftime')

 path.obs <- '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h-r1440x721cds/$var$_$date$.nc'
 lonlat_temp_st$obs <- CST_Start(dataset = path.obs,
                                 var = 'tas',
                                 date = unique(format(dates, '
                                 ftime = startR::values(dates), 
                                 ftime_across = 'date',
                                 ftime_var = 'ftime',
                                 merge_across_dims = TRUE,
                                 split_multiselected_dims = TRUE,
                                 lat = startR::values(list(latmin, latmax)),
                                 lat_reorder = startR::Sort(decreasing = TRUE),
                                 lon = startR::values(list(lonmin, lonmax)),
                                 lon_reorder = startR::CircularSort(0, 360),
                                 synonims = list(lon = c('lon', 'longitude'),
                                                 lat = c('lat', 'latitude'),
                                                 ftime = c('ftime', 'time')),
                                 transform = startR::CDORemapper,
                                 transform_extra_cells = 2,
                                 transform_params = list(grid = 'r360x181',
                                                         method = 'conservative'),
                                 transform_vars = c('lat', 'lon'),
                                 return_vars = list(lon = NULL,
                                                    lat = NULL,
                                                    ftime = 'date'),
                                 retrieve = TRUE)

 library(lubridate)
 dates_exp <- lonlat_temp_st$exp$attrs$Dates
 lonlat_temp_st$exp$attrs$Dates <- floor_date(ymd_hms(dates_exp), unit = "months")
 dim(lonlat_temp_st$exp$attrs$Dates) <- dim(dates_exp)
 
 dates_obs <- lonlat_temp_st$obs$attrs$Dates
 lonlat_temp_st$obs$attrs$Dates <- floor_date(ymd_hms(dates_obs), unit = "months")
 dim(lonlat_temp_st$obs$attrs$Dates) <- dim(dates_obs)

Author(s)

Nicolau Manubens nicolau.manubens@bsc.es

Print method for s2dv_cube objects

Description

This is an S3 method of the generic 'print' for the class 's2dv_cube'. When you will call 'print' on an 's2dv_cube' object, this method will display the content of the object in a clear and informative way.

Usage

## S3 method for class 's2dv_cube'
print(x, ...)

Arguments

Details

The object will be displayed following 's2dv_cube' class conventions. The top-level elements are: 'Data', a multidimensional array containing the object's data; 'Dimensions', the dimensions of the array; 'Coordinates', the array coordinates that match its dimensions, explicit coordinates have an asterisk (*) at the beginning while index coordinates do not; and 'Attributes', which contains all the metadata of the object. For more information about the 's2dv_cube', see s2dv_cube() and as.s2dv_cube() functions.

Value

No return value, called for side effects. It prints a summary of an 's2dv_cube' object to the console.

Creation of a 's2dv_cube' object

Description

This function allows to create an 's2dv_cube' object by passing information through its parameters. This function will be needed if the data hasn't been loaded using CST_Start or has been transformed with other methods. An 's2dv_cube' object has many different components including metadata. This function will allow to create 's2dv_cube' objects even if not all elements are defined and for each expected missed parameter a warning message will be returned.

Usage

s2dv_cube(
  data,
  coords = NULL,
  varName = NULL,
  metadata = NULL,
  Datasets = NULL,
  Dates = NULL,
  when = NULL,
  source_files = NULL,
  ...
)

Arguments

Value

The function returns an object of class 's2dv_cube' with the following elements in the structure:

• 'data', array with named dimensions;

• 'dims', named vector of the data dimensions;

• 'coords', list of named vectors with the coordinates corresponding to the dimensions of the data parameter;

• 'attrs', named list with elements:

  • 'Dates', array with named temporal dimensions of class 'POSIXct' from time values in the data;

  • 'Variable', has the following components:

    • 'varName', with the short name of the loaded variable as specified in the parameter 'var';

    • 'metadata', named list of elements with variable metadata. They can be from coordinates variables (e.g. longitude) or main variables (e.g. 'var');

  • 'Datasets', character strings indicating the names of the dataset;

  • 'source_files', a vector of character strings with complete paths to all the found files involved in loading the data;

  • 'when', a time stamp of the date issued by the Start() or Load() call to obtain the data;

  • 'load_parameters', it contains the components used in the arguments to load the data from Start() or Load() functions.

Author(s)

Perez-Zanon Nuria, nuria.perez@bsc.es

See Also

Load and CST_Start

Examples

exp_original <- 1:100
dim(exp_original) <- c(lat = 2, time = 10, lon = 5)
exp1 <- s2dv_cube(data = exp_original)
class(exp1)
coords <- list(lon = seq(-10, 10, 5), lat = c(45, 50))
exp2 <- s2dv_cube(data = exp_original, coords = coords) 
class(exp2)
metadata <- list(tas = list(level = '2m'))
exp3 <- s2dv_cube(data = exp_original, coords = coords,
                 varName = 'tas', metadata = metadata)
class(exp3)
Dates = as.POSIXct(paste0(rep("01", 10), rep("01", 10), 1990:1999), format = "%d%m%Y")
dim(Dates) <- c(time = 10)
exp4 <- s2dv_cube(data = exp_original, coords = coords,
                 varName = 'tas', metadata = metadata,
                 Dates = Dates)  
class(exp4)
exp5 <- s2dv_cube(data = exp_original, coords = coords,
                 varName = 'tas', metadata = metadata,
                 Dates = Dates, when = "2019-10-23 19:15:29 CET")  
class(exp5)
exp6 <- s2dv_cube(data = exp_original, coords = coords,
                 varName = 'tas', metadata = metadata,
                 Dates = Dates,
                 when = "2019-10-23 19:15:29 CET", 
                 source_files = c("/path/to/file1.nc", "/path/to/file2.nc"))
class(exp6)
exp7 <- s2dv_cube(data = exp_original, coords = coords,
                 varName = 'tas', metadata = metadata,
                 Dates = Dates,
                 when = "2019-10-23 19:15:29 CET", 
                 source_files = c("/path/to/file1.nc", "/path/to/file2.nc"),
                 Datasets = list(
                   exp1 = list(InitializationsDates = list(Member_1 = "01011990", 
                                                           Members = "Member_1"))))  
class(exp7)
dim(exp_original) <- c(dataset = 1, member = 1, time = 10, lat = 2, lon = 5)
exp8 <- s2dv_cube(data = exp_original, coords = coords,
                 varName = 'tas', metadata = metadata,
                 Dates = Dates, original_dates = Dates)  
class(exp8)

AEMET Training Training method (pre-downscaling) based on analogs: synoptic situations and significant predictors.

Description

This function caracterizes the synoptic situations in a past period based on low resolution reanalysis data (e.g, ERAInterim 1.5º x 1.5º) and an observational high resolution (HR) dataset (AEMET 5 km gridded daily precipitation and maximum and minimum temperature) (Peral et al., 2017)). The method uses three domains:

• peninsular Spain and Balearic Islands domain (5 km resolution): HR domain

• synoptic domain (low resolution): it should be centered over Iberian Peninsula and cover enough extension to detect as much synoptic situations as possible.

• extended domain (low resolution): it is an extension of the synoptic domain. It is used for 'slp_ext' parameter (see 'slp_lon' and 'slp_lat' below).

Usage

training_analogs(
  pred,
  slp_ext,
  lon,
  lat,
  slp_lon,
  slp_lat,
  var,
  HR_path,
  tdates
)

Arguments

Value

A matrix list (e.g. restrain) as a result of characterize the past synoptic situations and the significant predictors needed to downscale seasonal forecast variables. For precipitation the output includes:

• 'um': u component of geostrophic wind in all period (numeric matrix with [time, gridpoint] dimensions).

• 'vm': v component of geostrophic wind in all period (numeric matrix with [time,gridpoint] dimensions).

• 'nger': number of synoptic situations (integer).

• 'gu92': u component of geostrophic wind for each synoptic situation (numeric matrix with [nger,gridpoint] dimensions).

• 'gv92': v component of geostrophic wind for each synoptic situation (numeric matrix with [nger, gridpoint] dimensions).

• 'gu52': u component of wind at 500 hPa for each synotic situation (numeric matrix with [nger, gridpoint] dimensions).

• 'gv52': v component of wind at 500 hPa for each synotic situation (numeric matrix with [nger, gridpoint] dimensions).

• 'neni': number of reference centers where predictors are calculated (integer).

• 'vdmin': minimum distances between each HR gridpoint and the four nearest synoptic gridpoints (numeric matrix with [nptos,4] dimensions) (nptos = number of HR gridpoints).

• 'vref': four nearest synoptic gridpoints to each HR gridpoint (integer matrix with [nptos, 4] dimensions).

• 'ccm': multiple correlation coeficients (numeric matrix with [nger, nptos] dimensions) indices:

  • 'lab_pred': numeric labels of selected predictors (integer matrix with [nger,nptos,11,1] dimensions).

  • 'cor_pred': partial correlation of selected predictors (numeric matrix with [nger,nptos,11,2] dimensions).

For maximum and minimum temperature the output includes:

• 'um': u component of geostrophic wind in all training period (numeric matrix with [time,gridpoint] dimensions).

• 'vm': v component of geostrophic wind in all training period (numeric matrix with [time,gridpoint] dimensions).

• 'insol': insolation in all training period (numeric vector with [time] dimension).

• 'neni': number of reference centers where predictors are calculated (integer).

• 'vdmin': minimum distances between each HR gridpoint and the four nearest synoptic gridpoints (numeric matrix with [nptos,4] dimensions) (nptos = number of HR gridpoints).

• 'vref': four nearest synoptic gridpoints to each HR gridpoint (integer matrix with [nptos,4] dimensions).

The output can directly use as argument to 'CST_AnalogsPredictors' function (e.g. resdowns <- CST_AnalogsPredictors(...,restrain)).

Author(s)

Marta Dominguez Alonso - AEMET, mdomingueza@aemet.es

Nuria Perez-Zanon - BSC, nuria.perez@bsc.es