Newer
Older
#'@title Computing two dinamical proxies of the attractor in s2dv_cube.
#'
#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it}
#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it}
#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it}
#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr}
#'
#'@description This function computes two dinamical proxies of the attractor:
#'The local dimension (d) and the inverse of the persistence (theta) for an
#'These two parameters will be used as a condition for the computation of dynamical
#'scores to measure predictability and to compute bias correction conditioned by
#'the dynamics with the function DynBiasCorrection
#'Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in
#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019).
#' The hammam effect or how a warm ocean enhances large scale atmospheric predictability.
#' Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 "
#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017).
#' Dynamical proxies of North Atlantic predictability and extremes.
#' Scientific Reports, 7-41278, 2017.
#'
#'@param data a s2dv_cube object with the data to create the attractor. Must be a matrix with the timesteps in nrow
#'and the grids in ncol(dat(time,grids)
#
#'@param quanti list of arbitrary length of secondary grids. Each secondary grid is to
#' be provided as a list of length 2 with longitudes and latitudes
#'
#'@param ncores The number of cores to use in parallel computation
#'
#'@return dim and theta
#'
#'@examples
#'# Example 1: Computing the attractor using simple s2dv data
#'attractor <- CST_ProxiesAttractor(data = lonlat_data$obs, quanti = 0.6)
#'@export
CST_ProxiesAttractor <- function(data, quanti, ncores = NULL){
if (!inherits(data, 's2dv_cube')) {
stop("Parameter 'data' must be of the class 's2dv_cube', ",
"as output by CSTools::CST_Load.")
}
if (is.null(quanti)) {
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
data$data <- ProxiesAttractor(data = data$data, quanti = quanti, ncores = ncores)
return(data)
}
#'@rdname ProxiesAttractor
#'
#'@title Computing two dinamical proxies of the attractor.
#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it}
#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it}
#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it}
#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr}
#'
#'@description This function computes two dinamical proxies of the attractor:
#'The local dimension (d) and the inverse of the persistence (theta).
#'These two parameters will be used as a condition for the computation of dynamical
#'scores to measure predictability and to compute bias correction conditioned by
#'the dynamics with the function DynBiasCorrection.
#'Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in:
#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019).
#' The hammam effect or how a warm ocean enhances large scale atmospheric predictability.
#' Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 "
#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017).
#' Dynamical proxies of North Atlantic predictability and extremes.
#' Scientific Reports, 7-41278, 2017.
#'
#'@param data a multidimensional array with named dimensions to create the attractor. It requires a temporal dimension named 'time' and spatial dimensions called 'lat' and 'lon', or 'latitude' and 'longitude' or 'grid'.
#'@param quanti list of arbitrary length of secondary grids. Each secondary grid is to
#' be provided as a list of length 2 with longitudes and latitudes.
#'@param iplot a logical indicating if to plot results.
#'@param ncores The number of cores to use in parallel computation.
#'
#'@param iplot FALSE by default. If TRUE the function returns a simple plot
#'with three pannels on the top left a plot for local dimension 'd', on the top
#'right a plot for the inverse of the persistence 'theta', on the bottom a plot
#'for the 'attractor' plotted with the two properties 'd' (x axis) and 'theta'
#'(y axis)
#'
#'@return dim and theta
#'@import multiApply
#'
#'@examples
#'# Example 1: Computing the attractor using simple data
#'# Creating an example of matrix data(time,grids):
#'mat <- array(rnorm(36 * 40), c(time = 36, grid = 40))
#'qm <- 0.90 # imposing a threshold
#'Attractor <- ProxiesAttractor(data = mat, quanti = qm)
#'
#'@export
ProxiesAttractor <- function(data, quanti, iplot = FALSE, ncores = NULL){
if (is.null(data)) {
stop("Parameter 'data' cannot be NULL.")
}
if (is.null(quanti)) {
stop("Parameter 'quanti' is mandatory")
}
if (!is.logical(iplot) || is.null(iplot)) {
stop("Parameter 'iplot' is required and needs to be TRUE or FALSE.")
}
if (any(names(dim(data)) %in% 'sdate')) {
if (any(names(dim(data)) %in% 'ftime')) {
data <- MergeDims(data, merge_dims = c('ftime', 'sdate'),
rename_dim = 'time')
}
}
if (!(any(names(dim(data)) %in% 'time'))){
stop("Parameter 'data' must have a temporal dimension named 'time'.")
if (any(names(dim(data)) %in% 'lat')) {
if (any(names(dim(data)) %in% 'lon')) {
data <- MergeDims(data, merge_dims = c('lon', 'lat'),
rename_dim = 'grid')
}
}
if (any(names(dim(data)) %in% 'latitude')) {
if (any(names(dim(data)) %in% 'longitude')) {
data <- MergeDims(data, merge_dims = c('longitude', 'latitude'),
rename_dim = 'grid')
}
}
if(!(any(names(dim(data)) %in% 'grid'))){
stop("Parameter 'data' must have a spatial dimension named 'grid'.")
}
attractor <- Apply(data, target_dims = c('time', 'grid'),
fun = .proxiesattractor,
quanti = quanti , iplot = FALSE, ncores = ncores)
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
attractor <- lapply(attractor,
FUN = function(x, dimname){
names(dim(x))[dimname] <- 'time'
return(x)},
dimname = which(names(dim(attractor[[1]])) == 'dim2'))
return(list(dim = attractor$dim, theta = attractor$theta))
}
.proxiesattractor <- function(data, quanti, iplot = FALSE) {
# expected dimensions data: time and grid
logdista <- Apply(data, target_dims = 'grid',
fun = function(x, y){
-log(colMeans((y - as.vector(x))^2))},
y = t(data))[[1]]
#Computation of theta
Theta <- function(logdista, quanti){
#Compute the thheshold corresponding to the quantile
thresh <- quantile(logdista, quanti, na.rm = TRUE)
logdista[which(logdista == 'Inf')] <- NaN
Li <- which(as.vector(logdista) > as.numeric(thresh))
#Length of each cluster
Ti <- diff(Li)
N <- length(Ti)
q <- 1 - quanti
Si <- Ti - 1
Nc <- length(which(Si > 0))
N <- length(Ti)
theta <- (sum(q * Si) + N + Nc - sqrt(((sum(q * Si) + N + Nc)^2) -
8 * Nc * sum(q * Si))) / (2 * sum(q * Si))
#Sort the exceedances
logdista <- sort(logdista)
#Find all the Peaks over Thresholds.
findidx <- which(as.vector(logdista) > as.numeric(thresh))
if(length(findidx) < 1) {
stop("Parameter 'quanti' is too high for the length of the data provided.")
}
logextr <- logdista[findidx[[1]]:(length(logdista) - 1)]
#The inverse of the dimension is just the average of the exceedances
dim <- 1 /mean(as.numeric(logextr) - as.numeric(thresh))
return(list(dim = dim, theta = theta))
}
names(dim(logdista)) <- c('dim1', 'dim2')
proxies <- Apply(data = list(logdista = logdista),
target_dims = list('dim1'), fun = Theta, quanti = quanti,
ncores = ncores)
if(iplot == TRUE) {
time = c(1:length(proxies$theta))
layout(matrix(c(1, 3, 2, 3), 2, 2))
plot(time, proxies$dim, xlab = 'time', ylab = 'd',
main = 'local dimension', type = 'l')
plot(time, proxies$theta, xlab = 'time', ylab = 'theta', main = 'theta')
plot(proxies$dim, proxies$theta, col = 'blue',
main = "Proxies of the Attractor",
xlab = "local dimension", ylab = "theta", lwd = 8, 'p')
}
return(list(dim = proxies$dim, theta = proxies$theta))