Package 's2dv'

Description

The Absolute Mean Bias Skill Score is based on the Absolute Mean Error (Wilks, 2011) between the ensemble mean forecast and the observations. It measures the accuracy of the forecast in comparison with a reference forecast to assess whether the forecast presents an improvement or a worsening with respect to that reference. The Mean Bias Skill Score ranges between minus infinite and 1. Positive values indicate that the forecast has higher skill than the reference forecast, while negative values indicate that it has a lower skill. Examples of reference forecasts are the climatological forecast (average of the observations), a previous model version, or another model. It is computed as AbsBiasSS = 1 - AbsBias_exp / AbsBias_ref. The statistical significance is obtained based on a Random Walk test at the confidence level specified (DelSole and Tippett, 2016). If there is more than one dataset, the result will be computed for each pair of exp and obs data.

Usage

AbsBiasSS(
  exp,
  obs,
  ref = NULL,
  time_dim = "sdate",
  memb_dim = NULL,
  dat_dim = NULL,
  na.rm = FALSE,
  sig_method.type = "two.sided.approx",
  alpha = 0.05,
  ncores = NULL
)

Arguments

Value

References

Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7 DelSole and Tippett, 2016; https://doi.org/10.1175/MWR-D-15-0218.1

Examples

exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50))
ref <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50))
obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, sdate = 50))
biasSS1 <- AbsBiasSS(exp = exp, obs = obs, ref = ref, memb_dim = 'member')
biasSS2 <- AbsBiasSS(exp = exp, obs = obs, ref = NULL, memb_dim = 'member')

Description

Calculate the spatial anomaly correlation coefficient (ACC) for the ensemble mean of each model and the corresponding references over a spatial domain. It can return a forecast time series if the data contain forest time dimension, and also the ACC mean over one dimension, e.g., start date dimension. The domain of interest can be specified by providing the list of longitudes/ latitudes of the data together with the corners of the domain: lonlatbox = c(lonmin, lonmax, latmin, latmax). The data will be adjusted to have a spatial mean of zero, then area weighting is applied. The formula is referenced from Wilks (2011; section 7.6.4; https://doi.org/10.1016/B978-0-12-385022-5.00008-7).

Usage

ACC(
  exp,
  obs,
  dat_dim = NULL,
  lat_dim = "lat",
  lon_dim = "lon",
  avg_dim = "sdate",
  memb_dim = "member",
  lat = NULL,
  lon = NULL,
  lonlatbox = NULL,
  alpha = 0.05,
  pval = TRUE,
  sign = FALSE,
  conf = TRUE,
  conftype = "parametric",
  ncores = NULL
)

Arguments

Value

A list containing the numeric arrays:

References

Joliffe and Stephenson (2012). Forecast Verification: A Practitioner's Guide in Atmospheric Science. Wiley-Blackwell.; Wilks (2011; section 7.6.4; https://doi.org/10.1016/B978-0-12-385022-5.00008-7).

Examples

sampleData$mod <- Season(sampleData$mod, monini = 11, moninf = 12, monsup = 2)
sampleData$obs <- Season(sampleData$obs, monini = 11, moninf = 12, monsup = 2) 
clim <- Clim(sampleData$mod, sampleData$obs)
ano_exp <- Ano(sampleData$mod, clim$clim_exp)
ano_obs <- Ano(sampleData$obs, clim$clim_obs)
acc <- ACC(ano_exp, ano_obs, lat = sampleData$lat, dat_dim = 'dataset')
acc_bootstrap <- ACC(ano_exp, ano_obs, conftype = 'bootstrap', 
                    lat = sampleData$lat, dat_dim = 'dataset')
# Combine acc results for PlotACC
res <- array(c(acc$conf.lower, acc$acc, acc$conf.upper, acc$p.val), 
            dim = c(dim(acc$acc), 4))
res_bootstrap <- array(c(acc$acc_conf.lower, acc$acc, acc$acc_conf.upper, acc$p.val),
                      dim = c(dim(acc$acc), 4))
 
PlotACC(res, startDates)
PlotACC(res_bootstrap, startDates)

Description

The Atlantic Multidecadal Variability (AMV), also known as Atlantic Multidecadal Oscillation (AMO), is a mode of natural variability of the sea surface temperatures (SST) over the North Atlantic Ocean on multi-decadal time scales. The AMV index is computed as the difference of weighted-averaged SST anomalies over the North Atlantic region (0ºN-60ºN, 280ºE-360ºE) and the weighted-averaged SST anomalies over 60ºS-60ºN, 0ºE-360ºE (Trenberth & Dennis, 2005; Doblas-Reyes et al., 2013). If different members and/or datasets are provided, the climatology (used to calculate the anomalies) is computed individually for all of them.

Usage

AMV(
  data,
  data_lats,
  data_lons,
  type,
  lat_dim = "lat",
  lon_dim = "lon",
  mask = NULL,
  monini = 11,
  fmonth_dim = "fmonth",
  sdate_dim = "sdate",
  indices_for_clim = NULL,
  year_dim = "year",
  month_dim = "month",
  na.rm = TRUE,
  ncores = NULL
)

Arguments

Value

A numerical array with the AMV index with the same dimensions as data except the lat_dim, lon_dim and fmonth_dim (month_dim) in case of decadal predictions (historical simulations or observations). In case of decadal predictions, a new dimension 'fyear' is added.

Examples

## Observations or reanalyses
obs <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12))
lat <- seq(-90, 90, 10)
lon <- seq(0, 360, 10)
index_obs <- AMV(data = obs, data_lats = lat, data_lons = lon, type = 'obs')

## Historical simulations
hist <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5))
lat <- seq(-90, 90, 10)
lon <- seq(0, 360, 10)
index_hist <- AMV(data = hist, data_lats = lat, data_lons = lon, type = 'hist')

## Decadal predictions
dcpp <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5))
lat <- seq(-90, 90, 10)
lon <- seq(0, 360, 10)
index_dcpp <- AMV(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1)

Description

Create animations of maps in an equi-rectangular or stereographic projection, showing the anomalies, the climatologies, the mean InterQuartile Range, Maximum-Mininum, Standard Deviation, Median Absolute Deviation, the trends, the RMSE, the correlation or the RMSSS, between modelled and observed data along the forecast time (lead-time) for all input experiments and input observational datasets.

Usage

AnimateMap(
  var,
  lon,
  lat,
  toptitle = rep("", 11),
  sizetit = 1,
  units = "",
  monini = 1,
  freq = 12,
  msk95lev = FALSE,
  brks = NULL,
  cols = NULL,
  filled.continents = FALSE,
  lonmin = 0,
  lonmax = 360,
  latmin = -90,
  latmax = 90,
  intlon = 20,
  intlat = 30,
  drawleg = TRUE,
  subsampleg = 1,
  colNA = "white",
  equi = TRUE,
  fileout = c("output1_animvsltime.gif", "output2_animvsltime.gif",
    "output3_animvsltime.gif"),
  ...
)

Arguments

Details

Examples of input:

1. Outputs from clim (exp, obs, memb = FALSE): (nmod, nltime, nlat, nlon) or (nobs, nltime, nlat, nlon)

2. Model output from load/ano/smoothing: (nmod, nmemb, sdate, nltime, nlat, nlon) then passed through spread(var, posdim = 2, narm = TRUE) & mean1dim(var, posdim = 3, narm = TRUE) or through trend(mean1dim(var, 2), posTR = 2): (nmod, 3, nltime, nlat, nlon) animates average along start dates of IQR/MaxMin/SD/MAD across members or trends of the ensemble-mean computed accross the start dates.

3. model and observed output from load/ano/smoothing: (nmod, nmemb, sdate, nltime, nlat, nlon) & (nobs, nmemb, sdate, nltime, nlat, nlon) then averaged along members mean1dim(var_exp/var_obs, posdim = 2): (nmod, sdate, nltime, nlat, nlon) (nobs, sdate, nltime, nlat, nlon) then passed through corr(exp, obs, posloop = 1, poscor = 2) or RMS(exp, obs, posloop = 1, posRMS = 2): (nmod, nobs, 3, nltime, nlat, nlon) animates correlations or RMS between each exp & each obs against leadtime.

Examples

# See ?Load for explanations on the first part of this example
 ## Not run: 
data_path <- system.file('sample_data', package = 's2dv')
expA <- list(name = 'experiment', path = file.path(data_path,
            'model/$EXP_NAME$/$STORE_FREQ$_mean/$VAR_NAME$_3hourly',
            '$VAR_NAME$_$START_DATE$.nc'))
obsX <- list(name = 'observation', path = file.path(data_path,
            '$OBS_NAME$/$STORE_FREQ$_mean/$VAR_NAME$',
            '$VAR_NAME$_$YEAR$$MONTH$.nc'))

# Now we are ready to use Load().
startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101')
sampleData <- Load('tos', list(expA), list(obsX), startDates,
                  output = 'lonlat', latmin = 27, latmax = 48, 
                  lonmin = -12, lonmax = 40)
 
## End(Not run)
 
clim <- Clim(sampleData$mod, sampleData$obs, memb = FALSE)
 ## Not run: 
AnimateMap(clim$clim_exp, sampleData$lon, sampleData$lat,
 toptitle = "climatology of decadal prediction", sizetit = 1, 
 units = "degree", brks = seq(270, 300, 3), monini = 11, freq = 12, 
 msk95lev = FALSE, filled.continents = TRUE, intlon = 10, intlat = 10,
 fileout = 'clim_dec.gif')
 
## End(Not run)
# More examples in s2dverification but are deleted for now

Description

This function computes anomalies from a multidimensional data array and a climatology array.

Usage

Ano(data, clim, ncores = NULL)

Arguments

Value

An array with same dimensions as parameter 'data'.

Examples

# Load sample data as in Load() example:
example(Load)
clim <- Clim(sampleData$mod, sampleData$obs)
ano_exp <- Ano(sampleData$mod, clim$clim_exp)
ano_obs <- Ano(sampleData$obs, clim$clim_obs)
## Not run: 
PlotAno(ano_exp, ano_obs, startDates, 
       toptitle = 'Anomaly', ytitle = c('K', 'K', 'K'), 
       legends = 'ERSST', biglab = FALSE)

## End(Not run)

Description

Compute the anomalies from the arrays of the experimental and observational data output by subtracting the climatologies computed with a leave-one-out cross validation technique and a per-pair method (Garcia-Serrano and Doblas-Reyes, CD, 2012). Per-pair climatology means that only the start dates covered by the whole experiments/observational datasets will be used. In other words, the startdates which do not all have values along 'dat_dim' dimension of both the 'exp' and 'obs' are excluded when computing the climatologies.

Usage

Ano_CrossValid(
  exp,
  obs,
  time_dim = "sdate",
  dat_dim = c("dataset", "member"),
  memb_dim = "member",
  memb = TRUE,
  ncores = NULL
)

Arguments

Value

A list of 2:

Examples

# Load sample data as in Load() example:
example(Load)
anomalies <- Ano_CrossValid(sampleData$mod, sampleData$obs)
## Not run: 
PlotAno(anomalies$exp, anomalies$obs, startDates, 
       toptitle = paste('anomalies'), ytitle = c('K', 'K', 'K'), 
       legends = 'ERSST', biglab = FALSE)

## End(Not run)

Description

The Mean Bias or Mean Error (Wilks, 2011) is defined as the mean difference between the ensemble mean forecast and the observations. It is a deterministic metric. Positive values indicate that the forecasts are on average too high and negative values indicate that the forecasts are on average too low. It also allows to compute the Absolute Mean Bias or bias without temporal mean. If there is more than one dataset, the result will be computed for each pair of exp and obs data.

Usage

Bias(
  exp,
  obs,
  time_dim = "sdate",
  memb_dim = NULL,
  dat_dim = NULL,
  na.rm = FALSE,
  absolute = FALSE,
  time_mean = TRUE,
  alpha = 0.05,
  ncores = NULL
)

Arguments

Value

A numerical array of bias with dimensions c(nexp, nobs, the rest dimensions of 'exp' except 'time_dim' (if time_mean = T) and 'memb_dim'). 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 alpha is specified, and absolute is FALSE, the result is a list with two elements: the bias as described above and the significance as a logical array with the same dimensions.

References

Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7

Examples

exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50))
obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, sdate = 50))
bias <- Bias(exp = exp, obs = obs, memb_dim = 'member')
bias2 <- Bias(exp = exp, obs = obs, memb_dim = 'member', alpha = 0.01)
abs_bias <- Bias(exp = exp, obs = obs, memb_dim = 'member', absolute = TRUE, alpha = NULL)

Description

Compute the Brier score (BS) and the components of its standard decompostion with the two within-bin components described in Stephenson et al., (2008). It also returns the bias-corrected decomposition of the BS (Ferro and Fricker, 2012). BSS has the climatology as the reference forecast.

Usage

BrierScore(
  exp,
  obs,
  thresholds = seq(0.1, 0.9, 0.1),
  time_dim = "sdate",
  dat_dim = NULL,
  memb_dim = NULL,
  ncores = NULL
)

Arguments

Value

A list that contains:

The data type and dimensions of the items depend on if the input 'exp' and 'obs' are:
(a) Vectors
(b) Arrays with 'dat_dim' specified
(c) Arrays with no 'dat_dim' specified
Items 'rel', 'res', 'unc', 'bs', 'bs_check_res', 'bss_res', 'gres', 'bs_check_gres', 'bss_gres', 'rel_bias_corrected', 'gres_bias_corrected', 'unc_bias_corrected', and 'bss_bias_corrected' are (a) a number (b) an array with dimensions c(nexp, nobs, all the rest dimensions in 'exp' and 'obs' except 'time_dim' and 'memb_dim') (c) an array with dimensions of 'exp' and 'obs' except 'time_dim' and 'memb_dim'
Items 'nk', 'fkbar', and 'okbar' are (a) a vector of length of bin number determined by 'threshold' (b) an array with dimensions c(nexp, nobs, no. of bins, all the rest dimensions in 'exp' and 'obs' except 'time_dim' and 'memb_dim') (c) an array with dimensions c(no. of bin, all the rest dimensions in 'exp' and 'obs' except 'time_dim' and 'memb_dim')

References

Wilks (2006) Statistical Methods in the Atmospheric Sciences.
Stephenson et al. (2008). Two extra components in the Brier score decomposition. Weather and Forecasting, 23: 752-757.
Ferro and Fricker (2012). A bias-corrected decomposition of the BS. Quarterly Journal of the Royal Meteorological Society, DOI: 10.1002/qj.1924.

Examples

# Inputs are vectors
exp <- runif(10)
obs <- round(exp)
x <- BrierScore(exp, obs)

# Inputs are arrays
example(Load)
clim <- Clim(sampleData$mod, sampleData$obs)
ano_exp <- Ano(sampleData$mod, clim$clim_exp)
ano_obs <- Ano(sampleData$obs, clim$clim_obs)
bins_ano_exp <- ProbBins(ano_exp, thr = c(1/3, 2/3))
bins_ano_obs <- ProbBins(ano_obs, thr = c(1/3, 2/3))
res <- BrierScore(bins_ano_exp, MeanDims(bins_ano_obs, 'member'), memb_dim = 'member')

Description

This function takes as inputs a multidimensional array (optional), a vector or matrix of longitudes, a vector or matrix of latitudes, a destination grid specification, and the name of a method to be used to interpolate (one of those available in the 'remap' utility in CDO). The interpolated array is returned (if provided) together with the new longitudes and latitudes.

CDORemap() permutes by default the dimensions of the input array (if needed), splits it in chunks (CDO can work with data arrays of up to 4 dimensions), generates a file with the data of each chunk, interpolates it with CDO, reads it back into R and merges it into a result array. If no input array is provided, the longitude and latitude vectors will be transformed only. If the array is already on the desired destination grid, no transformation is performed (this behvaiour works only for lonlat and gaussian grids).

Any metadata attached to the input data array, longitudes or latitudes will be preserved or accordingly modified.

Usage

CDORemap(
  data_array = NULL,
  lons,
  lats,
  grid,
  method,
  avoid_writes = TRUE,
  crop = TRUE,
  force_remap = FALSE,
  write_dir = tempdir(),
  print_sys_msg = FALSE,
  ncores = NULL
)

Arguments

Value

A list with the following components:

Examples

## Not run: 
# Interpolating only vectors of longitudes and latitudes
lon <- seq(0, 360 - 360/50, length.out = 50)
lat <- seq(-90, 90, length.out = 25)
tas2 <- CDORemap(NULL, lon, lat, 't170grid', 'bil', TRUE)

# Minimal array interpolation
tas <- array(1:50, dim = c(25, 50))
names(dim(tas)) <- c('lat', 'lon')
lon <- seq(0, 360 - 360/50, length.out = 50)
lat <- seq(-90, 90, length.out = 25)
tas2 <- CDORemap(tas, lon, lat, 't170grid', 'bil', TRUE)

# Metadata can be attached to the inputs. It will be preserved and 
# accordignly modified.
tas <- array(1:50, dim = c(25, 50))
names(dim(tas)) <- c('lat', 'lon')
lon <- seq(0, 360 - 360/50, length.out = 50)
metadata <- list(lon = list(units = 'degrees_east'))
attr(lon, 'variables') <- metadata
lat <- seq(-90, 90, length.out = 25)
metadata <- list(lat = list(units = 'degrees_north'))
attr(lat, 'variables') <- metadata
metadata <- list(tas = list(dim = list(lat = list(len = 25,
                                                 vals = lat),
                                      lon = list(len = 50,
                                                 vals = lon)
                                     )))
attr(tas, 'variables') <- metadata
tas2 <- CDORemap(tas, lon, lat, 't170grid', 'bil', TRUE)

# Arrays of any number of dimensions in any order can be provided.
num_lats <- 25
num_lons <- 50
tas <- array(1:(10*num_lats*10*num_lons*10), 
            dim = c(10, num_lats, 10, num_lons, 10))
names(dim(tas)) <- c('a', 'lat', 'b', 'lon', 'c')
lon <- seq(0, 360 - 360/num_lons, length.out = num_lons)
metadata <- list(lon = list(units = 'degrees_east'))
attr(lon, 'variables') <- metadata
lat <- seq(-90, 90, length.out = num_lats)
metadata <- list(lat = list(units = 'degrees_north'))
attr(lat, 'variables') <- metadata
metadata <- list(tas = list(dim = list(a = list(),
                                      lat = list(len = num_lats,
                                                 vals = lat),
                                      b = list(),
                                      lon = list(len = num_lons,
                                                 vals = lon),
                                      c = list()
                                     )))
attr(tas, 'variables') <- metadata
tas2 <- CDORemap(tas, lon, lat, 't17grid', 'bil', TRUE)
# The step of permutation can be avoided but more intermediate file writes
# will be performed.
tas2 <- CDORemap(tas, lon, lat, 't17grid', 'bil', FALSE)

# If the provided array has the longitude or latitude dimension in the 
# right-most position, the same number of file writes will be performed,
# even if avoid_wrties = FALSE.
num_lats <- 25
num_lons <- 50
tas <- array(1:(10*num_lats*10*num_lons*10), 
            dim = c(10, num_lats, 10, num_lons))
names(dim(tas)) <- c('a', 'lat', 'b', 'lon')
lon <- seq(0, 360 - 360/num_lons, length.out = num_lons)
metadata <- list(lon = list(units = 'degrees_east'))
attr(lon, 'variables') <- metadata
lat <- seq(-90, 90, length.out = num_lats)
metadata <- list(lat = list(units = 'degrees_north'))
attr(lat, 'variables') <- metadata
metadata <- list(tas = list(dim = list(a = list(),
                                      lat = list(len = num_lats,
                                                 vals = lat),
                                      b = list(),
                                      lon = list(len = num_lons,
                                                 vals = lon)
                                     )))
attr(tas, 'variables') <- metadata
tas2 <- CDORemap(tas, lon, lat, 't17grid', 'bil', TRUE)
tas2 <- CDORemap(tas, lon, lat, 't17grid', 'bil', FALSE)

# An example of an interpolation from and onto a rectangular regular grid
num_lats <- 25
num_lons <- 50
tas <- array(1:(1*num_lats*num_lons), dim = c(num_lats, num_lons))
names(dim(tas)) <- c('y', 'x')
lon <- array(seq(0, 360 - 360/num_lons, length.out = num_lons), 
            dim = c(num_lons, num_lats))
metadata <- list(lon = list(units = 'degrees_east'))
names(dim(lon)) <- c('x', 'y')
attr(lon, 'variables') <- metadata
lat <- t(array(seq(-90, 90, length.out = num_lats), 
              dim = c(num_lats, num_lons)))
metadata <- list(lat = list(units = 'degrees_north'))
names(dim(lat)) <- c('x', 'y')
attr(lat, 'variables') <- metadata
tas2 <- CDORemap(tas, lon, lat, 'r100x50', 'bil')

# An example of an interpolation from an irregular grid onto a gaussian grid
num_lats <- 25
num_lons <- 50
tas <- array(1:(10*num_lats*10*num_lons*10), 
            dim = c(10, num_lats, 10, num_lons))
names(dim(tas)) <- c('a', 'j', 'b', 'i')
lon <- array(seq(0, 360 - 360/num_lons, length.out = num_lons), 
            dim = c(num_lons, num_lats))
metadata <- list(lon = list(units = 'degrees_east'))
names(dim(lon)) <- c('i', 'j')
attr(lon, 'variables') <- metadata
lat <- t(array(seq(-90, 90, length.out = num_lats), 
        dim = c(num_lats, num_lons)))
metadata <- list(lat = list(units = 'degrees_north'))
names(dim(lat)) <- c('i', 'j')
attr(lat, 'variables') <- metadata
tas2 <- CDORemap(tas, lon, lat, 't17grid', 'bil')

# Again, the dimensions can be in any order
num_lats <- 25
num_lons <- 50
tas <- array(1:(10*num_lats*10*num_lons), 
            dim = c(10, num_lats, 10, num_lons))
names(dim(tas)) <- c('a', 'j', 'b', 'i')
lon <- array(seq(0, 360 - 360/num_lons, length.out = num_lons), 
            dim = c(num_lons, num_lats))
names(dim(lon)) <- c('i', 'j')
lat <- t(array(seq(-90, 90, length.out = num_lats), 
              dim = c(num_lats, num_lons)))
names(dim(lat)) <- c('i', 'j')
tas2 <- CDORemap(tas, lon, lat, 't17grid', 'bil')
tas2 <- CDORemap(tas, lon, lat, 't17grid', 'bil', FALSE)
# It is ossible to specify an external NetCDF file as target grid reference
tas2 <- CDORemap(tas, lon, lat, 'external_file.nc', 'bil')

## End(Not run)

Description

This function computes per-pair climatologies for the experimental and observational data using one of the following methods:

1. per-pair method (Garcia-Serrano and Doblas-Reyes, CD, 2012 https://doi.org/10.1007/s00382-012-1413-1)

2. Kharin method (Kharin et al, GRL, 2012 https://doi.org/10.1029/2012GL052647)

3. Fuckar method (Fuckar et al, GRL, 2014 https://doi.org/10.1002/2014GL060815)

Per-pair climatology means that only the startdates covered by the whole experiments/observational dataset will be used. In other words, the startdates which are not all available along 'dat_dim' dimension of both the 'exp' and 'obs' are excluded when computing the climatologies. Kharin method is the linear trend bias correction method, and Fuckar method is the initial condition bias correction method. The two methods both do the per-pair correction beforehand.

Usage

Clim(
  exp,
  obs,
  time_dim = "sdate",
  dat_dim = c("dataset", "member"),
  method = "clim",
  ftime_dim = "ftime",
  memb = TRUE,
  memb_dim = "member",
  na.rm = TRUE,
  ncores = NULL
)

Arguments

Value

A list of 2:

Examples

# Load sample data as in Load() example:
example(Load)
clim <- Clim(sampleData$mod, sampleData$obs)
clim2 <- Clim(sampleData$mod, sampleData$obs, method = 'kharin', memb = FALSE)
## Not run: 
PlotClim(clim$clim_exp, clim$clim_obs, 
        toptitle = paste('sea surface temperature climatologies'), 
        ytitle = 'K', monini = 11, listexp = c('CMIP5 IC3'), 
        listobs = c('ERSST'), biglab = FALSE)

## End(Not run)

Description

Generates a colorblind friendly color palette with color ranges useful in climate temperature variable plotting.

Usage

clim.palette(palette = "bluered")

clim.colors(n, palette = "bluered")

Arguments

Examples

lims <- seq(-1, 1, length.out = 21)

ColorBar(lims, color_fun = clim.palette('redyellow'))

cols <- clim.colors(20)
ColorBar(lims, cols)

Description

Compute cluster centers and their time series of occurrences, with the K-means clustering method using Euclidean distance, of an array of input data with any number of dimensions that at least contain time_dim. Specifically, it partitions the array along time axis in K groups or clusters in which each space vector/array belongs to (i.e., is a member of) the cluster with the nearest center or centroid. This function is a wrapper of kmeans() and relies on the NbClust package (Charrad et al., 2014 JSS) to determine the optimal number of clusters used for K-means clustering if it is not provided by users.

Usage

Cluster(
  data,
  weights = NULL,
  time_dim = "sdate",
  space_dim = NULL,
  nclusters = NULL,
  index = "sdindex",
  ncores = NULL
)

Arguments

Value

A list containing:

References

Wilks, 2011, Statistical Methods in the Atmospheric Sciences, 3rd ed., Elsevire, pp 676.

Examples

# Generating synthetic data
a1 <- array(dim = c(200, 4))
mean1 <- 0
sd1 <- 0.3 

c0 <- seq(1, 200)
c1 <- sort(sample(x = 1:200, size = sample(x = 50:150, size = 1), replace = FALSE))
x1 <- c(1, 1, 1, 1)
for (i1 in c1) {
 a1[i1, ] <- x1 + rnorm(4, mean = mean1, sd = sd1)
}

c1p5 <- c0[!(c0 %in% c1)]
c2 <- c1p5[seq(1, length(c1p5), 2)] 
x2 <- c(2, 2, 4, 4)
for (i2 in c2) {
 a1[i2, ] <- x2 + rnorm(4, mean = mean1, sd = sd1)
}

c3 <- c1p5[seq(2, length(c1p5), 2)]
x3 <- c(3, 3, 1, 1)
for (i3 in c3) {
 a1[i3, ] <- x3 + rnorm(4, mean = mean1, sd = sd1)
}

# Computing the clusters
names(dim(a1)) <- c('sdate', 'space')
res1 <- Cluster(data = a1, weights = array(1, dim = dim(a1)[2]), nclusters = 3)
res2 <- Cluster(data = a1, weights = array(1, dim = dim(a1)[2]))

Description

Generates a color bar to use as colouring function for map plots and optionally draws it (horizontally or vertically) to be added to map multipanels or plots. It is possible to draw triangles at the ends of the colour bar to represent values that go beyond the range of interest. A number of options is provided to adjust the colours and the position and size of the components. The drawn colour bar spans a whole figure region and is compatible with figure layouts.

The generated colour bar consists of a set of breaks that define the length(brks) - 1 intervals to classify each of the values in each of the grid cells of a two-dimensional field. The corresponding grid cell of a given value of the field will be coloured in function of the interval it belongs to.

The only mandatory parameters are 'var_limits' or 'brks' (in its second format, see below).

Usage

ColorBar(
  brks = NULL,
  cols = NULL,
  vertical = TRUE,
  subsampleg = NULL,
  bar_limits = NULL,
  var_limits = NULL,
  triangle_ends = NULL,
  col_inf = NULL,
  col_sup = NULL,
  color_fun = clim.palette(),
  plot = TRUE,
  draw_ticks = TRUE,
  draw_separators = FALSE,
  triangle_ends_scale = 1,
  extra_labels = NULL,
  title = NULL,
  title_scale = 1,
  label_scale = 1,
  tick_scale = 1,
  extra_margin = rep(0, 4),
  label_digits = 4,
  ...
)

Arguments

Value

Examples

cols <- c("dodgerblue4", "dodgerblue1", "forestgreen", "yellowgreen", "white",
         "white", "yellow", "orange", "red", "saddlebrown")
lims <- seq(-1, 1, 0.2)
ColorBar(lims, cols)

Description

Composite a multi-dimensional array which contains two spatial and one temporal dimensions, e.g., (lon, lat, time), according to the indices of mode/cluster occurrences in time. The p-value by t-test is also computed.

Usage

Composite(
  data,
  occ,
  time_dim = "time",
  space_dim = c("lon", "lat"),
  lag = 0,
  eno = FALSE,
  K = NULL,
  fileout = NULL,
  ncores = NULL
)

Arguments

Value

A list containing:

Examples

blank <- array(0, dim = c(20, 10, 30))
x1 <- blank
t1 <- blank
f1 <- blank

for (i in 1:20) {
 x1[i, , ] <- i
}

for (i in 1:30) {
 t1[, , i] <- i
}

# This is 2D propagating sin wave example, where we use f1(lon, lat, time) 
# wave field. Compositing (like using stroboscopicc light) at different
# time steps can lead to modification or cancelation of wave pattern.

for (i in 1:20) {
 for (j in 1:30) {
   f1[i, , j] <- 3 * sin(2 * pi * x1[i, , j] / 5. - 2 * pi * t1[i, , j] / 6.)
 }
}
names(dim(f1)) <- c('lon', 'lat', 'time')
occ <- rep(0, 30)
occ[c(2, 5, 8, 11, 14, 17, 20, 23)] <- 1
res <- Composite(data = f1, occ = occ)
filled.contour(res$composite[, , 1])

occ <- rep(0, 30)
occ[c(3, 9, 15, 21)] <- 1
res <- Composite(data = f1, occ = occ)
filled.contour(res$composite[, , 1])

# Example with one missing composite in occ:
data <- 1:(4 * 5 * 6)
dim(data) <- c(lon = 4, lat = 5, case = 6)
occ <- c(1, 1, 2, 2, 3, 3) 
res <- Composite(data, occ, time_dim = 'case',  K = 4)

Description

Given a pair of dataset name and variable name, this function determines applies all the matching entries found in the corresponding configuration table to work out the dataset main path, file path, actual name of variable inside NetCDF files, ...

Usage

ConfigApplyMatchingEntries(
  configuration,
  var,
  exp = NULL,
  obs = NULL,
  show_entries = FALSE,
  show_result = TRUE
)

Arguments

Value

A list with the information resulting of applying the matching entries is returned.

See Also

ConfigApplyMatchingEntries, ConfigEditDefinition, ConfigEditEntry, ConfigFileOpen, ConfigShowSimilarEntries, ConfigShowTable

Examples

# Create an empty configuration file
config_file <- paste0(tempdir(), "/example.conf")
s2dv::ConfigFileCreate(config_file, confirm = FALSE)
# Open it into a configuration object
configuration <- ConfigFileOpen(config_file)
# Add an entry at the bottom of 4th level of file-per-startdate experiments 
# table which will associate the experiment "ExampleExperiment2" and variable 
# "ExampleVariable" to some information about its location.
configuration <- ConfigAddEntry(configuration, "experiments", 
                "last", "ExampleExperiment2", "ExampleVariable", 
                "/path/to/ExampleExperiment2/", 
                "ExampleVariable/ExampleVariable_$START_DATE$.nc")
# Edit entry to generalize for any variable. Changing variable needs .
configuration <- ConfigEditEntry(configuration, "experiments", 1, 
                var_name = ".*", 
                file_path = "$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc")
# Now apply matching entries for variable and experiment name and show the 
# result
match_info <- ConfigApplyMatchingEntries(configuration, 'tas', 
             exp = c('ExampleExperiment2'), show_result = TRUE)

Description

These functions help in adding, modifying or removing variable definitions in a configuration object obtained with ConfigFileOpen or ConfigFileCreate. ConfigEditDefinition() will add the definition if not existing.

Usage

ConfigEditDefinition(configuration, name, value, confirm = TRUE)

ConfigRemoveDefinition(configuration, name)

Arguments

Value

A modified configuration object is returned.

See Also

[ConfigApplyMatchingEntries()], [ConfigEditDefinition()], [ConfigEditEntry()], [ConfigFileOpen()], [ConfigShowSimilarEntries()], [ConfigShowTable()].

Examples

# Create an empty configuration file
config_file <- paste0(tempdir(), "/example.conf")
ConfigFileCreate(config_file, confirm = FALSE)
# Open it into a configuration object
configuration <- ConfigFileOpen(config_file)
# Add an entry at the bottom of 4th level of file-per-startdate experiments 
# table which will associate the experiment "ExampleExperiment2" and variable 
# "ExampleVariable" to some information about its location.
configuration <- ConfigAddEntry(configuration, "experiments", 
                "last", "ExampleExperiment2", "ExampleVariable", 
                "/path/to/ExampleExperiment2/", 
                "ExampleVariable/ExampleVariable_$START_DATE$.nc")
# Edit entry to generalize for any variable. Changing variable needs .
configuration <- ConfigEditEntry(configuration, "experiments", 1, 
                var_name = ".*", 
                file_path = "$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc")
# Now apply matching entries for variable and experiment name and show the 
# result
match_info <- ConfigApplyMatchingEntries(configuration, 'tas', 
             exp = c('ExampleExperiment2'), show_result = TRUE)

Description

ConfigAddEntry(), ConfigEditEntry() and ConfigRemoveEntry() are functions to manage entries in a configuration object created with ConfigFileOpen().
Before adding an entry, make sure the defaults don't do already what you want (ConfigShowDefinitions(), ConfigShowTable()).
Before adding an entry, make sure it doesn't override and spoil what other entries do (ConfigShowTable(), ConfigFileOpen()).
Before adding an entry, make sure there aren't other entries that already do what you want (ConfigShowSimilarEntries()).

Usage

ConfigEditEntry(
  configuration,
  dataset_type,
  position,
  dataset_name = NULL,
  var_name = NULL,
  main_path = NULL,
  file_path = NULL,
  nc_var_name = NULL,
  suffix = NULL,
  varmin = NULL,
  varmax = NULL
)

ConfigAddEntry(
  configuration,
  dataset_type,
  position = "last",
  dataset_name = ".*",
  var_name = ".*",
  main_path = "*",
  file_path = "*",
  nc_var_name = "*",
  suffix = "*",
  varmin = "*",
  varmax = "*"
)

ConfigRemoveEntry(
  configuration,
  dataset_type,
  dataset_name = NULL,
  var_name = NULL,
  position = NULL
)

Arguments

Value

The function returns an accordingly modified configuration object. To apply the changes in the configuration file it must be saved using ConfigFileSave().

See Also

ConfigApplyMatchingEntries, ConfigEditDefinition, ConfigEditEntry, ConfigFileOpen, ConfigShowSimilarEntries, ConfigShowTable

Examples

# Create an empty configuration file
config_file <- paste0(tempdir(), "/example.conf")
ConfigFileCreate(config_file, confirm = FALSE)
# Open it into a configuration object
configuration <- ConfigFileOpen(config_file)
# Add an entry at the bottom of 4th level of file-per-startdate experiments 
# table which will associate the experiment "ExampleExperiment" and variable 
# "ExampleVariable" to some information about its location.
configuration <- ConfigAddEntry(configuration, "experiments", 
                "last", "ExampleExperiment", "ExampleVariable", 
                "/path/to/ExampleExperiment/", 
                "ExampleVariable/ExampleVariable_$START_DATE$.nc")
# Add another entry
configuration <- ConfigAddEntry(configuration, "experiments",
                "last", "ExampleExperiment2", "ExampleVariable", 
                "/path/to/ExampleExperiment2/",
                "ExampleVariable/ExampleVariable_$START_DATE$.nc")
# Edit second entry to generalize for any variable. Changing variable needs .
configuration <- ConfigEditEntry(configuration, "experiments", 2, 
                var_name = ".*",
                file_path = "$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc")
# Remove first entry
configuration <- ConfigRemoveEntry(configuration, "experiments",
                "ExampleExperiment", "ExampleVariable")
# Show results
ConfigShowTable(configuration, "experiments")
# Save the configuration
ConfigFileSave(configuration, config_file, confirm = FALSE)

Description

These functions help in creating, opening and saving configuration files.

Usage

ConfigFileOpen(file_path, silent = FALSE, stop = FALSE)

ConfigFileCreate(file_path, confirm = TRUE)

ConfigFileSave(configuration, file_path, confirm = TRUE)

Arguments

Details

ConfigFileOpen() loads all the data contained in the configuration file specified as parameter 'file_path'. Returns a configuration object with the variables needed for the configuration file mechanism to work. This function is called from inside the Load() function to load the configuration file specified in 'configfile'.

ConfigFileCreate() creates an empty configuration file and saves it to the specified path. It may be opened later with ConfigFileOpen() to be edited. Some default values are set when creating a file with this function, you can check these with ConfigShowDefinitions().

ConfigFileSave() saves a configuration object into a file, which may then be used from Load().

Two examples of configuration files can be found inside the 'inst/config/' folder in the package:

• BSC.conf: configuration file used at BSC-CNS. Contains location data on several datasets and variables.

• template.conf: very simple configuration file intended to be used as pattern when starting from scratch.

How the configuration file works:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~
It contains one list and two tables.
Each of these have a header that starts with '!!'. These are key lines and should not be removed or reordered.
Lines starting with '#' and blank lines will be ignored. The list should contains variable definitions and default value definitions.
The first table contains information about experiments.
The third table contains information about observations.
Each table entry is a list of comma-separated elements.
The two first are part of a key that is associated to a value formed by the other elements.
The key elements are a dataset identifier and a variable name.
The value elements are the dataset main path, dataset file path, the variable name inside the .nc file, a default suffix (explained below) and a minimum and maximum vaues beyond which loaded data is deactivated.
Given a dataset name and a variable name, a full path is obtained concatenating the main path and the file path.
Also the nc variable name, the suffixes and the limit values are obtained.
Any of the elements in the keys can contain regular expressions[1] that will cause matching for sets of dataset names or variable names.
The dataset path and file path can contain shell globbing expressions[2] that will cause matching for sets of paths when fetching the file in the full path.
The full path can point to an OPeNDAP URL.
Any of the elements in the value can contain variables that will be replaced to an associated string.
Variables can be defined only in the list at the top of the file.
The pattern of a variable definition is
VARIABLE_NAME = VARIABLE_VALUE
and can be accessed from within the table values or from within the variable values as
$VARIABLE_NAME$
For example:
FILE_NAME = tos.nc
!!table of experiments
ecmwf, tos, /path/to/dataset/, $FILE_NAME$
There are some reserved variables that will offer information about the store frequency, the current startdate Load() is fetching, etc:
$VAR_NAME$, $START_DATE$, $STORE_FREQ$, $MEMBER_NUMBER$
for experiments only: $EXP_NAME$
for observations only: $OBS_NAME$, $YEAR$, $MONTH$, $DAY$
Additionally, from an element in an entry value you can access the other elements of the entry as:
$EXP_MAIN_PATH$, $EXP_FILE_PATH$,
$VAR_NAME$, $SUFFIX$, $VAR_MIN$, $VAR_MAX$

The variable $SUFFIX$ is useful because it can be used to take part in the main or file path. For example: '/path/to$SUFFIX$/dataset/'.
It will be replaced by the value in the column that corresponds to the suffix unless the user specifies a different suffix via the parameter 'suffixexp' or 'suffixobs'.
This way the user is able to load two variables with the same name in the same dataset but with slight modifications, with a suffix anywhere in the path to the data that advices of this slight modification.

The entries in a table will be grouped in 4 levels of specificity:

1. General entries:
   - the key dataset name and variable name are both a regular expression matching any sequence of characters (.*) that will cause matching for any pair of dataset and variable names
   Example: .*, .*, /dataset/main/path/, file/path, nc_var_name, suffix, var_min, var_max

2. Dataset entries:
   - the key variable name matches any sequence of characters
   Example: ecmwf, .*, /dataset/main/path/, file/path, nc_var_name, suffix, var_min, var_max

3. Variable entries:
   - the key dataset name matches any sequence of characters
   Example: .*, tos, /dataset/main/path/, file/path, nc_var_name, suffix, var_min, var_max

4. Specific entries:
   - both key values are specified
   Example: ecmwf, tos, /dataset/main/path/, file/path, nc_var_name, suffix, var_min, var_max

Given a pair of dataset name and variable name for which we want to know the full path, all the rules that match will be applied from more general to more specific.
If there is more than one entry per group that match a given key pair, these will be applied in the order of appearance in the configuration file (top to bottom).

An asterisk (*) in any value element will be interpreted as 'leave it as is or take the default value if yet not defined'.
The default values are defined in the following reserved variables:
$DEFAULT_EXP_MAIN_PATH$, $DEFAULT_EXP_FILE_PATH$, $DEFAULT_NC_VAR_NAME$, $DEFAULT_OBS_MAIN_PATH$, $DEFAULT_OBS_FILE_PATH$, $DEFAULT_SUFFIX$, $DEFAULT_VAR_MIN$, $DEFAULT_VAR_MAX$,
$DEFAULT_DIM_NAME_LATITUDES$, $DEFAULT_DIM_NAME_LONGITUDES$,
$DEFAULT_DIM_NAME_MEMBERS$

Trailing asterisks in an entry are not mandatory. For example
ecmwf, .*, /dataset/main/path/, *, *, *, *, *
will have the same effect as
ecmwf, .*, /dataset/main/path/

A double quote only (") in any key or value element will be interpreted as 'fill in with the same value as the entry above'.

Value

ConfigFileOpen() returns a configuration object with all the information for the configuration file mechanism to work.
ConfigFileSave() returns TRUE if the file has been saved and FALSE otherwise.
ConfigFileCreate() returns nothing.

References

[1] https://stat.ethz.ch/R-manual/R-devel/library/base/html/regex.html
[2] https://tldp.org/LDP/abs/html/globbingref.html

See Also

ConfigApplyMatchingEntries, ConfigEditDefinition, ConfigEditEntry, ConfigFileOpen, ConfigShowSimilarEntries, ConfigShowTable

Examples

# Create an empty configuration file
config_file <- paste0(tempdir(), "/example.conf")
ConfigFileCreate(config_file, confirm = FALSE)
# Open it into a configuration object
configuration <- ConfigFileOpen(config_file)
# Add an entry at the bottom of 4th level of file-per-startdate experiments 
# table which will associate the experiment "ExampleExperiment2" and variable 
# "ExampleVariable" to some information about its location.
configuration <- ConfigAddEntry(configuration, "experiments", 
                "last", "ExampleExperiment2", "ExampleVariable", 
                "/path/to/ExampleExperiment2/", 
                "ExampleVariable/ExampleVariable_$START_DATE$.nc")
# Edit entry to generalize for any variable. Changing variable needs .
configuration <- ConfigEditEntry(configuration, "experiments", 1, 
                var_name = ".*", 
                file_path = "$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc")
# Now apply matching entries for variable and experiment name and show the 
# result
match_info <- ConfigApplyMatchingEntries(configuration, 'tas', 
             exp = c('ExampleExperiment2'), show_result = TRUE)
# Finally save the configuration file.
ConfigFileSave(configuration, config_file, confirm = FALSE)

Description

These functions help in finding similar entries in tables of supported datasets by comparing all entries with some given information.
This is useful when dealing with complex configuration files and not sure if already support certain variables or datasets.
At least one field must be provided in ConfigShowSimilarEntries(). Other fields can be unspecified and won't be taken into account. If more than one field is provided, sameness is avreaged over all provided fields and entries are sorted from higher average to lower.

Usage

ConfigShowSimilarEntries(
  configuration,
  dataset_name = NULL,
  var_name = NULL,
  main_path = NULL,
  file_path = NULL,
  nc_var_name = NULL,
  suffix = NULL,
  varmin = NULL,
  varmax = NULL,
  n_results = 10
)

Arguments

Details

Sameness is calculated with string distances as specified by Simon White in [1].

Value

These functions return information about the found matches.

References

[1] Simon White, string seamness: http://www.catalysoft.com/articles/StrikeAMatch.html

See Also

ConfigApplyMatchingEntries, ConfigEditDefinition, ConfigEditEntry, ConfigFileOpen, ConfigShowSimilarEntries, ConfigShowTable

Examples

# Create an empty configuration file
config_file <- paste0(tempdir(), "/example.conf")
ConfigFileCreate(config_file, confirm = FALSE)
# Open it into a configuration object
configuration <- ConfigFileOpen(config_file)
# Add an entry at the bottom of 4th level of file-per-startdate experiments 
# table which will associate the experiment "ExampleExperiment2" and variable 
# "ExampleVariable" to some information about its location.
configuration <- ConfigAddEntry(configuration, "experiments", "last", 
                "ExampleExperiment2", "ExampleVariable", 
                "/path/to/ExampleExperiment2/", 
                "ExampleVariable/ExampleVariable_$START_DATE$.nc")
# Edit entry to generalize for any variable. Changing variable needs .
configuration <- ConfigEditEntry(configuration, "experiments", 1, 
                var_name = "Var.*", 
                file_path = "$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc")
# Look for similar entries
ConfigShowSimilarEntries(configuration, dataset_name = "Exper", 
                        var_name = "Vari")

Description

These functions show the tables of supported datasets and definitions in a configuration object obtained via ConfigFileCreate() or ConfigFileOpen().

Usage

ConfigShowTable(configuration, dataset_type, line_numbers = NULL)

ConfigShowDefinitions(configuration)

Arguments

Value

These functions return nothing.

See Also

[ConfigApplyMatchingEntries()], [ConfigEditDefinition()], [ConfigEditEntry()], [ConfigFileOpen()], [ConfigShowSimilarEntries()], [ConfigShowTable()].

Examples

# Create an empty configuration file
config_file <- paste0(tempdir(), "/example.conf")
ConfigFileCreate(config_file, confirm = FALSE)
# Open it into a configuration object
configuration <- ConfigFileOpen(config_file)
# Add an entry at the bottom of 4th level of file-per-startdate experiments 
# table which will associate the experiment "ExampleExperiment2" and variable 
# "ExampleVariable" to some information about its location.
configuration <- ConfigAddEntry(configuration, "experiments", "last", 
                "ExampleExperiment2", "ExampleVariable", 
                "/path/to/ExampleExperiment2/", 
                "ExampleVariable/ExampleVariable_$START_DATE$.nc")
# Edit entry to generalize for any variable. Changing variable needs .
configuration <- ConfigEditEntry(configuration, "experiments", 1, 
                var_name = ".*", 
                file_path = "$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc")
# Show tables, lists and definitions
ConfigShowTable(configuration, 'experiments')
ConfigShowDefinitions(configuration)

Description

Compute the linear trend for a time series by least square fitting together with the associated error interval for both the observational and model data. The 95% confidence interval and detrended observational and model data are also provided.
The function doesn't do the ensemble mean, so if the input data have the member dimension, ensemble mean needs to be computed beforehand.

Usage

Consist_Trend(
  exp,
  obs,
  dat_dim = "dataset",
  time_dim = "sdate",
  interval = 1,
  ncores = NULL
)

Arguments

Value

A list containing:

Examples

#'# Load sample data as in Load() example:
example(Load)
clim <- Clim(sampleData$mod, sampleData$obs)
ano_exp <- Ano(sampleData$mod, clim$clim_exp)
ano_obs <- Ano(sampleData$obs, clim$clim_obs)
runmean_months <- 12
smooth_ano_exp <- Smoothing(ano_exp, runmeanlen = runmean_months)
smooth_ano_obs <- Smoothing(ano_obs, runmeanlen = runmean_months)
dim_to_mean <- 'member'  # average along members
years_between_startdates <- 5
trend <- Consist_Trend(MeanDims(smooth_ano_exp, dim_to_mean, na.rm = TRUE), 
                      MeanDims(smooth_ano_obs, dim_to_mean, na.rm = TRUE), 
                      interval = years_between_startdates)
#Bind data for plotting
trend_bind <- abind::abind(trend$conf.lower[2, , ], trend$trend[2, , ], 
                          trend$conf.upper[2, , ], trend$trend[1, , ], along = 0)
trend_bind <- Reorder(trend_bind, c(2, 1, 3))

PlotVsLTime(trend_bind, toptitle = "trend", ytitle = "K/(5 years)", 
           monini = 11, limits = c(-0.8, 0.8), listexp = c('CMIP5 IC3'), 
           listobs = c('ERSST'), biglab = FALSE, hlines = c(0)) 
PlotAno(InsertDim(trend$detrended_exp, 2, 1), InsertDim(trend$detrended_obs, 2, 1), 
       startDates, "Detrended tos anomalies", ytitle = 'K', 
       legends = 'ERSST', biglab = FALSE)

Description

Calculate the correlation coefficient (Pearson, Kendall or Spearman) for an array of forecast and an array of observation. The correlations are computed along 'time_dim' that usually refers to the start date dimension. If 'comp_dim' is given, the correlations are computed only if obs along comp_dim dimension are complete between limits[1] and limits[2], i.e., there is no NA between limits[1] and limits[2]. This option can be activated if the user wants to account only for the forecasts which the corresponding observations are available at all leadtimes.
The confidence interval is computed by the Fisher transformation and the significance level relies on an one-sided student-T distribution.
The function can calculate ensemble mean before correlation by 'memb_dim' specified and 'memb = F'. If ensemble mean is not calculated, correlation will be calculated for each member. If there is only one dataset for exp and obs, you can simply use cor() to compute the correlation.

Usage

Corr(
  exp,
  obs,
  time_dim = "sdate",
  dat_dim = NULL,
  comp_dim = NULL,
  limits = NULL,
  method = "pearson",
  memb_dim = NULL,
  memb = TRUE,
  pval = TRUE,
  conf = TRUE,
  sign = FALSE,
  alpha = 0.05,
  ncores = NULL
)

Arguments

Value

A list containing the numeric arrays with dimension:
c(nexp, nobs, exp_memb, obs_memb, all other dimensions of exp except time_dim and memb_dim).
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. exp_memb is the number of member in experiment (i.e., 'memb_dim' in exp) and obs_memb is the number of member in observation (i.e., 'memb_dim' in obs). If memb = F, exp_memb and obs_memb are omitted.

Examples

# Case 1: Load sample data as in Load() example: 
example(Load) 
clim <- Clim(sampleData$mod, sampleData$obs) 
ano_exp <- Ano(sampleData$mod, clim$clim_exp) 
ano_obs <- Ano(sampleData$obs, clim$clim_obs) 
runmean_months <- 12 

# Smooth along lead-times   
smooth_ano_exp <- Smoothing(ano_exp, runmeanlen = runmean_months) 
smooth_ano_obs <- Smoothing(ano_obs, runmeanlen = runmean_months) 
required_complete_row <- 3  # Discard start dates which contain any NA lead-times 
leadtimes_per_startdate <- 60 
corr <- Corr(MeanDims(smooth_ano_exp, 'member'),              
            MeanDims(smooth_ano_obs, 'member'),              
            comp_dim = 'ftime', dat_dim = 'dataset', 
            limits = c(ceiling((runmean_months + 1) / 2),                         
            leadtimes_per_startdate - floor(runmean_months / 2))) 

# Case 2: Keep member dimension
corr <- Corr(smooth_ano_exp, smooth_ano_obs, memb_dim = 'member', dat_dim = 'dataset')
# ensemble mean
corr <- Corr(smooth_ano_exp, smooth_ano_obs, memb_dim = 'member', memb = FALSE,
            dat_dim = 'dataset')

Description

The Continuous Ranked Probability Score (CRPS; Wilks, 2011) is the continuous version of the Ranked Probability Score (RPS; Wilks, 2011). It is a skill metric to evaluate the full distribution of probabilistic forecasts. It has a negative orientation (i.e., the higher-quality forecast the smaller CRPS) and it rewards the forecast that has probability concentration around the observed value. In case of a deterministic forecast, the CRPS is reduced to the mean absolute error. It has the same units as the data. The function is based on enscrps_cpp from SpecsVerification. If there is more than one dataset, CRPS will be computed for each pair of exp and obs data.

Usage

CRPS(
  exp,
  obs,
  time_dim = "sdate",
  memb_dim = "member",
  dat_dim = NULL,
  Fair = FALSE,
  return_mean = TRUE,
  ncores = NULL
)

Arguments

Value

A numerical array of CRPS with dimensions c(nexp, nobs, the rest dimensions of 'exp' except 'time_dim' and 'memb_dim' dimensions). 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.

References

Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7

Examples

exp <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
obs <- array(rnorm(1000), dim = c(lat = 3, lon = 2, sdate = 50))
res <- CRPS(exp = exp, obs = obs)

Description

The Continuous Ranked Probability Skill Score (CRPSS; Wilks, 2011) is the skill score based on the Continuous Ranked Probability Score (CRPS; Wilks, 2011). It can be used to assess whether a forecast presents an improvement or worsening with respect to a reference forecast. The CRPSS ranges between minus infinite and 1. If the CRPSS is positive, it indicates that the forecast has higher skill than the reference forecast, while a negative value means that it has a lower skill. Examples of reference forecasts are the climatological forecast, persistence, a previous model version, or another model. It is computed as 'CRPSS = 1 - CRPS_exp / CRPS_ref. The statistical significance is obtained based on a Random Walk test at the specified confidence level (DelSole and Tippett, 2016).

Usage

CRPSS(
  exp,
  obs,
  ref = NULL,
  time_dim = "sdate",
  memb_dim = "member",
  dat_dim = NULL,
  Fair = FALSE,
  clim.cross.val = TRUE,
  sig_method.type = "two.sided.approx",
  alpha = 0.05,
  N.eff = NA,
  ncores = NULL
)

Arguments

Value

References

Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7 DelSole and Tippett, 2016; https://doi.org/10.1175/MWR-D-15-0218.1

Examples

exp <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
obs <- array(rnorm(1000), dim = c(lat = 3, lon = 2, sdate = 50))
ref <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
res <- CRPSS(exp = exp, obs = obs) ## climatology as reference forecast
res <- CRPSS(exp = exp, obs = obs, ref = ref) ## ref as reference forecast

Description

Compute the correlation difference between two deterministic forecasts. Positive values of the correlation difference indicate that the forecast is more skillful than the reference forecast, while negative values mean that the reference forecast is more skillful. The statistical significance of the correlation differences is computed with a one-sided or two-sided test for equality of dependent correlation coefficients (Steiger, 1980; Siegert et al., 2017) using effective degrees of freedom to account for the autocorrelation of the time series (Zwiers and von Storch, 1995).

Usage

DiffCorr(
  exp,
  obs,
  ref,
  N.eff = NA,
  time_dim = "sdate",
  memb_dim = NULL,
  method = "pearson",
  alpha = 0.05,
  handle.na = "return.na",
  test.type = "two-sided",
  pval = TRUE,
  sign = FALSE,
  ncores = NULL
)

Arguments

Value

A list with:

References

Steiger, 1980; https://content.apa.org/doi/10.1037/0033-2909.87.2.245 Siegert et al., 2017; https://doi.org/10.1175/MWR-D-16-0037.1 Zwiers and von Storch, 1995; https://doi.org/10.1175/1520-0442(1995)008<0336:TSCIAI>2.0.CO;2

Examples

exp <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
obs <- array(rnorm(1000), dim = c(lat = 3, lon = 2, sdate = 50))
ref <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
res_two.sided_sign <- DiffCorr(exp, obs, ref, memb_dim = 'member', 
                               test.type = 'two-sided', alpha = 0.05)
res_one.sided_pval <- DiffCorr(exp, obs, ref, memb_dim = 'member', 
                               test.type = 'one-sided', alpha = NULL)

Description

Compute the number of effective samples along one dimension of an array. This effective number of independent observations can be used in statistical/inference tests.
The calculation is based on eno function from Caio Coelho from rclim.txt.

Usage

Eno(data, time_dim = "sdate", na.action = na.pass, ncores = NULL)

Arguments

Value

An array with the same dimension as parameter 'data' except the time_dim dimension, which is removed after the computation. The array indicates the number of effective sample along time_dim.

Examples

set.seed(1)
data <- array(rnorm(800), dim = c(dataset = 1, member = 2, sdate = 4, 
                                 ftime = 4, lat = 10, lon = 10))
na <- floor(runif(40, min = 1, max = 800))
data[na] <- NA
res <- Eno(data)

Description

Perform an area-weighted EOF analysis using single value decomposition (SVD) based on a covariance matrix or a correlation matrix if parameter 'corr' is set to TRUE.

Usage

EOF(
  ano,
  lat,
  lon,
  time_dim = "sdate",
  space_dim = c("lat", "lon"),
  neofs = 15,
  corr = FALSE,
  ncores = NULL
)

Arguments

Value

A list containing:

See Also

ProjectField, NAO, PlotBoxWhisker

Examples

# This example computes the EOFs along forecast horizons and plots the one 
# that explains the greatest amount of variability. The example data has low  
# resolution so the result may not be explanatory, but it displays how to 
# use this function.

ano <- Ano_CrossValid(sampleData$mod, sampleData$obs)
tmp <- MeanDims(ano$exp, c('dataset', 'member'))
ano <- tmp[1, , ,]
names(dim(ano)) <- names(dim(tmp))[-2]
eof <- EOF(ano, sampleData$lat, sampleData$lon)
## Not run: 
PlotEquiMap(eof$EOFs[1, , ], sampleData$lon, sampleData$lat)

## End(Not run)

Description

Calculate the four main teleconnection indices in European Atlantic Ocean region: North Atlantic oscillation (NAO), East Atlantic Pattern (EA), East Atlantic/Western Russia (EAWR), and Scandinavian pattern (SCA). The function REOF() is used for the calculation, and the first four modes are returned.

Usage

EuroAtlanticTC(
  ano,
  lat,
  lon,
  ntrunc = 30,
  time_dim = "sdate",
  space_dim = c("lat", "lon"),
  corr = FALSE,
  ncores = NULL
)

Arguments

Value

A list containing:

See Also

REOF NAO

Examples

# Use synthetic data
set.seed(1)
dat <- array(rnorm(800), dim = c(dat = 2, sdate = 5, lat = 8, lon = 15))
lat <- seq(10, 90, length.out = 8)
lon <- seq(-100, 70, length.out = 15)
res <- EuroAtlanticTC(dat, lat = lat, lon = lon)

Description

Filter out the selected frequency from a time series. The filtering is performed by dichotomy, seeking for a frequency around the parameter 'freq' and the phase that maximizes the signal to subtract from the time series. The maximization of the signal to subtract relies on a minimization of the mean square differences between the time series ('data') and the cosine of the specified frequency and phase.

Usage

Filter(data, freq, time_dim = "ftime", ncores = NULL)

Arguments

Value

A numeric vector or array of the filtered data with the dimensions the same as 'data'.

Examples

# Load sample data as in Load() example:
example(Load)
ensmod <- MeanDims(sampleData$mod, 2)
spectrum <- Spectrum(ensmod)

for (jsdate in 1:dim(spectrum)['sdate']) {
 for (jlen in 1:dim(spectrum)['ftime']) {
   if (spectrum[jlen, 2, 1, jsdate] > spectrum[jlen, 3, 1, jsdate]) {
     ensmod[1, jsdate, ] <- Filter(ensmod[1, jsdate, ], spectrum[jlen, 1, 1, jsdate])
   }
 }
}
 
PlotAno(InsertDim(ensmod, 2, 1), sdates = startDates)

Description

Compute probabilistic forecasts from an ensemble based on the relative thresholds, or the probabilistic observations (i.e., which probabilistic category was observed). A reference period can be specified to calculate the absolute thresholds between each probabilistic category. The absolute thresholds can be computed in cross-validation mode. If data is an ensemble, the probabilities are calculated as the percentage of members that fall into each category. For observations (or forecast without member dimension), 1 means that the event happened, while 0 indicates that the event did not happen. Weighted probabilities can be computed if the weights are provided for each ensemble member and time step. The absolute thresholds can also be provided directly for probabilities calculation.

Usage

GetProbs(
  data,
  time_dim = "sdate",
  memb_dim = "member",
  indices_for_quantiles = NULL,
  prob_thresholds = c(1/3, 2/3),
  abs_thresholds = NULL,
  bin_dim_abs = "bin",
  weights = NULL,
  cross.val = FALSE,
  ncores = NULL
)

Arguments

Value

A numerical array of probabilities with dimensions c(bin_dim_abs, the rest dimensions of 'data' except 'memb_dim'). 'bin' dimension has the length of probabilistic categories, i.e., length(prob_thresholds) + 1.

Examples

data <- array(rnorm(2000), dim = c(ensemble = 25, sdate = 20, time = 4))
res <- GetProbs(data = data, time_dim = 'sdate', memb_dim = 'ensemble', 
               indices_for_quantiles = 4:17)

# abs_thresholds is provided
abs_thr1 <- c(-0.2, 0.3)
abs_thr2 <- array(c(-0.2, 0.3) + rnorm(40) * 0.1, dim = c(cat = 2, sdate = 20))
res1 <- GetProbs(data = data, time_dim = 'sdate', memb_dim = 'ensemble',
                prob_thresholds = NULL, abs_thresholds = abs_thr1)
res2 <- GetProbs(data = data, time_dim = 'sdate', memb_dim = 'ensemble',
                prob_thresholds = NULL, abs_thresholds = abs_thr2, bin_dim_abs = 'cat')

Description

The Global Mean Surface Temperature (GMST) anomalies are computed as the weighted-averaged surface air temperature anomalies over land and sea surface temperature anomalies over the ocean. If different members and/or datasets are provided, the climatology (used to calculate the anomalies) is computed individually for all of them.

Usage

GMST(
  data_tas,
  data_tos,
  data_lats,
  data_lons,
  mask_sea_land,
  sea_value,
  type,
  mask = NULL,
  lat_dim = "lat",
  lon_dim = "lon",
  monini = 11,
  fmonth_dim = "fmonth",
  sdate_dim = "sdate",
  indices_for_clim = NULL,
  year_dim = "year",
  month_dim = "month",
  na.rm = TRUE,
  ncores = NULL
)

Arguments

Value

A numerical array with the GMST anomalies with the same dimensions as data_tas except the lat_dim, lon_dim and fmonth_dim (month_dim) in case of decadal predictions (historical simulations or observations). In case of decadal predictions, a new dimension 'fyear' is added.

Examples

## Observations or reanalyses
obs_tas <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12))
obs_tos <- array(2:101, dim = c(year = 5, lat = 19, lon = 37, month = 12))
mask_sea_land <- array(c(1,0,1), dim = c(lat = 19, lon = 37))
sea_value <- 1
lat <- seq(-90, 90, 10)
lon <- seq(0, 360, 10)
index_obs <- GMST(data_tas = obs_tas, data_tos = obs_tos, data_lats = lat,
                  data_lons = lon, type = 'obs', 
                  mask_sea_land = mask_sea_land, sea_value = sea_value)

## Historical simulations
hist_tas <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5))
hist_tos <- array(2:101, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5))
mask_sea_land <- array(c(1,0,1), dim = c(lat = 19, lon = 37))
sea_value <- 1
lat <- seq(-90, 90, 10)
lon <- seq(0, 360, 10)
index_hist <- GMST(data_tas = hist_tas, data_tos = hist_tos, data_lats = lat, 
                   data_lons = lon, type = 'hist', mask_sea_land = mask_sea_land, 
                   sea_value = sea_value)

## Decadal predictions
dcpp_tas <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5))
dcpp_tos <- array(2:101, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5))
mask_sea_land <- array(c(1,0,1), dim = c(lat = 19, lon = 37))
sea_value <- 1
lat <- seq(-90, 90, 10)
lon <- seq(0, 360, 10)
index_dcpp <- GMST(data_tas = dcpp_tas, data_tos = dcpp_tos, data_lats = lat, 
                   data_lons = lon, type = 'dcpp', monini = 1, mask_sea_land = mask_sea_land,
                   sea_value = sea_value)

Description

The Global Surface Air Temperature (GSAT) anomalies are computed as the weighted-averaged surface air temperature anomalies over the global region. If different members and/or datasets are provided, the climatology (used to calculate the anomalies) is computed individually for all of them.

Usage

GSAT(
  data,
  data_lats,
  data_lons,
  type,
  lat_dim = "lat",
  lon_dim = "lon",
  mask = NULL,
  monini = 11,
  fmonth_dim = "fmonth",
  sdate_dim = "sdate",
  indices_for_clim = NULL,
  year_dim = "year",
  month_dim = "month",
  na.rm = TRUE,
  ncores = NULL
)

Arguments

Value

A numerical array with the GSAT anomalies with the same dimensions as data except the lat_dim, lon_dim and fmonth_dim (month_dim) in case of decadal predictions (historical simulations or observations). In case of decadal predictions, a new dimension 'fyear' is added.

Examples

## Observations or reanalyses
obs <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12))
lat <- seq(-90, 90, 10)
lon <- seq(0, 360, 10)
index_obs <- GSAT(data = obs, data_lats = lat, data_lons = lon, type = 'obs')

## Historical simulations
hist <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5))
lat <- seq(-90, 90, 10)
lon <- seq(0, 360, 10)
index_hist <- GSAT(data = hist, data_lats = lat, data_lons = lon, type = 'hist')

## Decadal predictions
dcpp <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5))
lat <- seq(-90, 90, 10)
lon <- seq(0, 360, 10)
index_dcpp <- GSAT(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1)

Description

Reorganize a long run (historical typically) with only one start date into chunks corresponding to a set of start dates. The time frequency of the data should be monthly.

Usage

Histo2Hindcast(
  data,
  sdatesin,
  sdatesout,
  nleadtimesout,
  sdate_dim = "sdate",
  ftime_dim = "ftime",
  ncores = NULL
)

Arguments

Value

A numeric array with the same dimensions as data, except the length of sdate_dim is 'sdatesout' and the length of ftime_dim is nleadtimesout.

Examples

sdates_out <- c('19901101', '19911101', '19921101', '19931101', '19941101')
leadtimes_per_startdate <- 12	
exp_data <- Histo2Hindcast(sampleData$mod, startDates, 
                          sdates_out, leadtimes_per_startdate)
obs_data <- Histo2Hindcast(sampleData$obs, startDates, 
                          sdates_out, leadtimes_per_startdate)
 ## Not run: 
exp_data <- Reorder(exp_data, c(3, 4, 1, 2))
obs_data <- Reorder(obs_data, c(3, 4, 1, 2))
PlotAno(exp_data, obs_data, sdates_out, 
       toptitle = paste('Anomalies reorganized into shorter chunks'), 
       ytitle = 'K', fileout = NULL)
 
## End(Not run)

Description

Insert an extra dimension into an array at position 'posdim' with length 'lendim'. The array repeats along the new dimension.

Usage

InsertDim(data, posdim, lendim, name = NULL)

Arguments

Value

An array as parameter 'data' but with the added named dimension.

Examples

a <- array(rnorm(15), dim = c(a = 3, b = 1, c = 5, d = 1))
res <- InsertDim(InsertDim(a, posdim = 2, lendim = 1, name = 'e'), 4, c(f = 2))
dim(res)