Package 'CSTools'

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, [email protected] for PROSNOW version

Lauriane Batté, [email protected] 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'))

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, [email protected]

Maria M. Chaves-Montero, [email protected]

Veronica Torralba, [email protected]

Nuria Perez-Zanon [email protected]

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. [email protected]

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)

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, [email protected]

Nicolau Manubens, [email protected]

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)

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, [email protected]

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)

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, [email protected]

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)

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, [email protected]

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)

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, [email protected]

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)

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, [email protected]

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

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, [email protected]

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)

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-one-out cross validation) calibration are possible.

Usage

Calibration(
  exp,
  obs,
  exp_cor = NULL,
  cal.method = "mse_min",
  eval.method = "leave-one-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
)

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, [email protected]

Bert Van Schaeybroeck, [email protected]

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)

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, ...)

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, [email protected]

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

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, [email protected] for PROSNOW version

Lauriane Batté, [email protected] 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)

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é, [email protected]

Paola Marson, [email protected]

Gildas Dayon, [email protected]

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'))

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, [email protected]

Maria M. Chaves-Montero, [email protected]

Veronica Torralba, [email protected]

Nuria Perez-Zanon [email protected]

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. [email protected]

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)

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, [email protected]

Nuria Perez-Zanon - BSC, [email protected]

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, [email protected]

Pena Jesus, [email protected]

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)

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, [email protected]

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)

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, [email protected]

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)

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-one-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
)

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, [email protected]

Bert Van Schaeybroeck, [email protected]

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")

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-one-out",
  amt.cat = 3,
  ...
)

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, [email protected]

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")

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, [email protected]

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)

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, [email protected]

Maria M. Chaves-Montero, [email protected]

Veronica Torralba, [email protected]

Davide Faranda, [email protected]

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)

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 = F
)

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, [email protected]

Ignazio Giuntoli - ISAC-CNR, [email protected]

Danila Volpi - ISAC-CNR, [email protected]

Paolo Davini - ISAC-CNR, [email protected]

Jost von Hardenberg - ISAC-CNR, [email protected]

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"))

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, [email protected]

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

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, [email protected]

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)

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

Author(s)

Nuria Perez-Zanon, [email protected]

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)

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, [email protected]

Paolo Davini - ISAC-CNR, [email protected]

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)

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, [email protected]

Perez-Zanon Nuria, [email protected]

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)

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, [email protected]

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))

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 ([email protected]) used in

Usage

CST_ProxiesAttractor(data, quanti, ncores = NULL)

Arguments

Value

dim and theta

Author(s)

Carmen Alvarez-Castro, [email protected]

Maria M. Chaves-Montero, [email protected]

Veronica Torralba, [email protected]

Davide Faranda, [email protected]

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)

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,
  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, [email protected]

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)

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, [email protected]

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')

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