Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
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
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
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
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
#'Compute the North Atlantic Oscillation (NAO) Index
#'
#'Compute the North Atlantic Oscillation (NAO) index based on the leading EOF
#'of the sea level pressure (SLP) anomalies over the north Atlantic region
#'(20N-80N, 80W-40E). The PCs are obtained by projecting the forecast and
#'observed anomalies onto the observed EOF pattern (Pobs) or the forecast
#'anomalies onto the EOF pattern of the other years of the forecast (Pmod).
#'By default (ftime_avg = 2:4) NAO() computes the NAO index for 1-month
#'lead seasonal forecasts that can be plotted with PlotBoxWhisker(). It returns
#'cross-validated PCs of the NAO index for forecast (exp) and observations
#'(obs) based on the leading EOF pattern.
#'
#'@param exp A named numeric array of North Atlantic SLP (20N-80N, 80W-40E)
#' forecast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with
#' dimensions 'time_dim', 'memb_dim', 'ftime_dim', and 'space_dim' at least.
#' If only NAO of observational data needs to be computed, this parameter can
#' be left to NULL. The default value is NULL.
#'@param obs A named numeric array of North Atlantic SLP (20N-80N, 80W-40E)
#' observed anomalies from \code{Ano()} or \code{Ano_CrossValid()} with
#' dimensions 'time_dim', 'memb_dim', 'ftime_dim', and 'space_dim' at least.
#' If only NAO of experimental data needs to be computed, this parameter can
#' be left to NULL. The default value is NULL.
#'@param lat A vector of the latitudes of 'exp' and 'obs'.
#'@param lon A vector of the longitudes of 'exp' and 'obs'.
#'@param time_dim A character string indicating the name of the time dimension
#' of 'exp' and 'obs'. The default value is 'sdate'.
#'@param memb_dim A character string indicating the name of the member
#' dimension of 'exp' and 'obs'. The default value is 'member'.
#'@param space_dim A vector of two character strings. The first is the dimension
#' name of latitude of 'ano' and the second is the dimension name of longitude
#' of 'ano'. The default value is c('lat', 'lon').
#'@param ftime_dim A character string indicating the name of the forecast time
#' dimension of 'exp' and 'obs'. The default value is 'ftime'.
#'@param ftime_avg A numeric vector of the forecast time steps to average
#' across the target period. The default value is 2:4, i.e., from 2nd to 4th
#' forecast time steps.
#'@param obsproj A logical value indicating whether to compute the NAO index by
#' projecting the forecast anomalies onto the leading EOF of observational
#' reference (TRUE) or compute the NAO by first computing the leading
#' EOF of the forecast anomalies (in cross-validation mode, i.e. leaving the
#' year you are evaluating out), and then projecting forecast anomalies onto
#' this EOF (FALSE). The default value is TRUE.
#'@param ncores An integer indicating the number of cores to use for parallel
#' computation. The default value is NULL.
#'
#'@return
#'A list which contains:
#'\item{exp}{
#' A numeric array of forecast NAO index in verification format with the same
#' dimensions as 'exp' except space_dim and ftime_dim.
#' }
#'\item{obs}{
#' A numeric array of observed NAO index in verification format with the same
#' dimensions as 'obs' except space_dim and ftime_dim.
#'}
#'
#'@references
#'Doblas-Reyes, F.J., Pavan, V. and Stephenson, D. (2003). The skill of
#' multi-model seasonal forecasts of the wintertime North Atlantic
#' Oscillation. Climate Dynamics, 21, 501-514.
#' DOI: 10.1007/s00382-003-0350-4
#'
#'@examples
#' \dontshow{
#'startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101')
#'sampleData <- s2dv:::.LoadSampleData('tos', c('experiment'),
#' c('observation'), startDates,
#' leadtimemin = 1,
#' leadtimemax = 4,
#' output = 'lonlat',
#' latmin = 27, latmax = 48,
#' lonmin = -12, lonmax = 40)
#'# No example data is available over NAO region, so in this example we will
#'# tweak the available data. In a real use case, one can Load() the data over
#'# the NAO region directly.
#'sampleData$lon[] <- c(40, 280, 340)
#'sampleData$lat[] <- c(20, 80)
#' }
#'
#'# Now ready to compute the EOFs and project on, for example, the first
#'# variability mode.
#'ano <- Ano_CrossValid(sampleData$mod, sampleData$obs)
#'# Note that computing the NAO over the region for which there is available
#'# example data is not the full NAO area: NAO() will raise a warning.
#'nao <- NAO(ano$exp, ano$obs, sampleData$lat, sampleData$lon)
#'# Finally plot the NAO index
#' \dontrun{
#'nao$exp <- Reorder(nao$exp, c(2, 1))
#'nao$obs <- Reorder(nao$obs, c(2, 1))
#'PlotBoxWhisker(nao$exp, nao$obs, "NAO index, DJF", "NAO index (PC1) TOS",
#' monini = 12, yearini = 1985, freq = 1, "Exp. A", "Obs. X")
#' }
#'
#'@import multiApply
#'@importFrom ClimProjDiags Subset
#'@export
NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
memb_dim = 'member', space_dim = c('lat', 'lon'),
ftime_dim = 'ftime', ftime_avg = 2:4,
obsproj = TRUE, ncores = NULL) {
# Check inputs
## exp and obs (1)
if (is.null(obs) & is.null(exp)) {
stop("Parameter 'exp' and 'obs' cannot both be NULL.")
}
if (!is.null(exp)) {
if (!is.numeric(exp)) {
stop("Parameter 'exp' must be a numeric array.")
}
if (is.null(dim(exp))) {
stop(paste0("Parameter 'exp' and must have at least dimensions ",
"time_dim, memb_dim, space_dim, and ftime_dim."))
}
if(any(is.null(names(dim(exp)))) | any(nchar(names(dim(exp))) == 0)) {
stop("Parameter 'exp' must have dimension names.")
}
}
if (!is.null(obs)) {
if (!is.numeric(obs)) {
stop("Parameter 'obs' must be a numeric array.")
}
if (is.null(dim(obs))) {
stop(paste0("Parameter 'obs' and must have at least dimensions ",
"time_dim, memb_dim, space_dim, and ftime_dim."))
}
if(any(is.null(names(dim(obs)))) | any(nchar(names(dim(obs))) == 0)) {
stop("Parameter 'obs' must have dimension names.")
}
}
## time_dim
if (!is.character(time_dim) | length(time_dim) > 1) {
stop("Parameter 'time_dim' must be a character string.")
}
if (!is.null(exp)) {
if (!time_dim %in% names(dim(exp))) {
stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.")
}
}
if (!is.null(obs)) {
if (!time_dim %in% names(dim(obs))) {
stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.")
}
}
## memb_dim
if (!is.character(memb_dim) | length(memb_dim) > 1) {
stop("Parameter 'memb_dim' must be a character string.")
}
if (!is.null(exp)) {
if (!memb_dim %in% names(dim(exp))) {
stop("Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension.")
}
}
if (!is.null(obs)) {
if (!memb_dim %in% names(dim(obs))) {
stop("Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension.")
}
}
## space_dim
if (!is.character(space_dim) | length(space_dim) != 2) {
stop("Parameter 'space_dim' must be a character vector of 2.")
}
if (!is.null(exp)) {
if (any(!space_dim %in% names(dim(exp)))) {
stop("Parameter 'space_dim' is not found in 'exp' or 'obs' dimension.")
}
}
if (!is.null(obs)) {
if (any(!space_dim %in% names(dim(obs)))) {
stop("Parameter 'space_dim' is not found in 'exp' or 'obs' dimension.")
}
}
## ftime_dim
if (!is.character(ftime_dim) | length(ftime_dim) > 1) {
stop("Parameter 'ftime_dim' must be a character string.")
}
if (!is.null(exp)) {
if (!ftime_dim %in% names(dim(exp))) {
stop("Parameter 'ftime_dim' is not found in 'exp' or 'obs' dimension.")
}
}
if (!is.null(obs)) {
if (!ftime_dim %in% names(dim(obs))) {
stop("Parameter 'ftime_dim' is not found in 'exp' or 'obs' dimension.")
}
}
## exp and obs (2)
if (!is.null(exp) & !is.null(obs)) {
name_exp <- sort(names(dim(exp)))
name_obs <- sort(names(dim(obs)))
name_exp <- name_exp[-which(name_exp == memb_dim)]
name_obs <- name_obs[-which(name_obs == memb_dim)]
if(!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
stop(paste0("Parameter 'exp' and 'obs' must have the same length of ",
"all dimensions except 'memb_dim'."))
}
}
## ftime_avg
if (!is.vector(ftime_avg) | !is.integer(ftime_avg)) {
stop("Parameter 'ftime_avg' must be an integer vector.")
}
if (!is.null(exp)) {
if (max(ftime_avg) > dim(exp)[ftime_dim] | min(ftime_avg) < 1) {
stop("Parameter 'ftime_avg' must be within the range of ftime_dim length.")
}
} else {
if (max(ftime_avg) > dim(obs)[ftime_dim] | min(ftime_avg) < 1) {
stop("Parameter 'ftime_avg' must be within the range of ftime_dim length.")
}
}
## sdate >= 2
if (!is.null(exp)) {
if (dim(exp)[time_dim] < 2) {
stop("The length of time_dim must be at least 2.")
}
} else {
if (dim(obs)[time_dim] < 2) {
stop("The length of time_dim must be at least 2.")
}
}
## lat and lon
if (!is.null(exp)) {
if (!is.numeric(lat) | length(lat) != dim(exp)[space_dim[1]]) {
stop(paste0("Parameter 'lat' must be a numeric vector with the same ",
"length as the latitude dimension of 'exp' and 'obs'."))
}
if (!is.numeric(lon) | length(lon) != dim(exp)[space_dim[2]]) {
stop(paste0("Parameter 'lon' must be a numeric vector with the same ",
"length as the longitude dimension of 'exp' and 'obs'."))
}
} else {
if (!is.numeric(lat) | length(lat) != dim(obs)[space_dim[1]]) {
stop(paste0("Parameter 'lat' must be a numeric vector with the same ",
"length as the latitude dimension of 'exp' and 'obs'."))
}
if (!is.numeric(lon) | length(lon) != dim(obs)[space_dim[2]]) {
stop(paste0("Parameter 'lon' must be a numeric vector with the same ",
"length as the longitude dimension of 'exp' and 'obs'."))
}
}
stop_needed <- FALSE
if (tail(lat, 1) < 70 | tail(lat, 1) > 90 |
head(lat, 1) > 30 | head(lat, 1) < 10) {
stop_needed <- TRUE
}
#NOTE: different from s2dverification
# lon is not used in the calculation actually. EOF only uses lat to do the
# weight. So we just need to ensure the data is in this region, regardless
# the order.
if (any(lon < 0)) { #[-180, 180]
if (!(min(lon) > -90 & min(lon) < -70 & max(lon) < 50 & max(lon) > 30)) {
stop_needed <- TRUE
}
} else { #[0, 360]
if (any(lon >= 50 & lon <= 270)) {
stop_needed <- TRUE
} else {
lon_E <- lon[which(lon < 50)]
lon_W <- lon[-which(lon < 50)]
if (max(lon_E) < 30 | min(lon_W) > 290) {
stop_needed <- TRUE
}
}
}
if (stop_needed) {
stop(paste0("The typical domain used to compute the NAO is 20N-80N, ",
"80W-40E. 'lat' or 'lon' is out of range."))
}
## obsproj
if (!is.logical(obsproj) | length(obsproj) > 1) {
stop("Parameter 'obsproj' must be either TRUE or FALSE.")
}
if (obsproj) {
if (is.null(obs)) {
stop("Parameter 'obsproj' set to TRUE but no 'obs' provided.")
}
if (is.null(exp)) {
.warning("parameter 'obsproj' set to TRUE but no 'exp' provided.")
}
}
## ncores
if (!is.null(ncores)) {
if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
length(ncores) > 1) {
stop("Parameter 'ncores' must be a positive integer.")
}
}
#-------- Average ftime -----------
if (!is.null(exp)) {
exp_sub <- ClimProjDiags::Subset(exp, ftime_dim, ftime_avg, drop = FALSE)
exp <- MeanDims(exp_sub, ftime_dim, na.rm = TRUE)
## Cross-validated PCs. Fabian. This should be extended to
## nmod and nlt by simple loops. Virginie
}
if (!is.null(obs)) {
obs_sub <- ClimProjDiags::Subset(obs, ftime_dim, ftime_avg, drop = FALSE)
obs <- MeanDims(obs_sub, ftime_dim, na.rm = TRUE)
}
if (!is.null(exp) & !is.null(obs)) {
res <- Apply(list(exp, obs),
target_dims = list(exp = c(memb_dim, time_dim, space_dim),
obs = c(memb_dim, time_dim, space_dim)),
fun = .NAO,
obsproj = obsproj, lat = lat, lon = lon,
ncores = ncores)
} else if (!is.null(exp)) {
res <- Apply(list(exp = exp),
target_dims = list(exp = c(memb_dim, time_dim, space_dim)),
fun = .NAO,
obsproj = obsproj, lat = lat, lon = lon, obs = NULL,
ncores = ncores)
} else if (!is.null(obs)) {
res <- Apply(list(obs = obs),
target_dims = list(obs = c(memb_dim, time_dim, space_dim)),
fun = .NAO,
obsproj = obsproj, lat = lat, lon = lon, exp = NULL,
ncores = ncores)
}
return(res)
}
.NAO <- function(exp = NULL, obs = NULL, lat, lon,
obsproj = TRUE, ncores = NULL) {
# exp: [memb_exp, sdate, lat, lon]
# obs: [memb_obs, sdate, lat, lon]
if (!is.null(exp)) {
ntime <- dim(exp)[2]
nlat <- dim(exp)[3]
nlon <- dim(exp)[4]
nmemb_exp <- dim(exp)[1]
nmemb_obs <- dim(obs)[1]
} else {
ntime <- dim(obs)[2]
nlat <- dim(obs)[3]
nlon <- dim(obs)[4]
nmemb_obs <- dim(obs)[1]
}
if (!is.null(obs)) NAOO.ver <- array(NA, dim = c(ntime, nmemb_obs))
if (!is.null(exp)) NAOF.ver <- array(NA, dim = c(ntime, nmemb_exp))
for (tt in 1:ntime) { #sdate
if (!is.null(obs)) {
## Observed EOF excluding one forecast start year.
obs_sub <- ClimProjDiags::Subset(obs, 2, c(1:ntime)[-tt], drop = FALSE)
obs_EOF <- EOF(obs_sub, lat = lat, lon = lon, time_dim = names(ntime),
space_dim = c(names(nlat), names(nlon)), neofs = 1)
## Correct polarity of pattern.
#NOTE: different from s2dverification
# dim(obs_EOF$EOFs): [mode, lat, lon, member]
for (imemb in 1:nmemb_obs) {
if (0 < mean(obs_EOF$EOFs[1, which.min(abs(lat - 65)), , ], na.rm = T)) {
obs_EOF$EOFs[1, , , imemb] <- obs_EOF$EOFs[1, , , imemb] * (-1)
}
}
# obs_EOF$PCs <- obs_EOF$PCs * sign # not used
## Project observed anomalies.
PF <- ProjectField(obs, eof = obs_EOF, time_dim = names(ntime),
space_dim = c(names(nlat), names(nlon)), mode = 1)
NAOO.ver[tt, ] <- PF[tt, ]
## Keep PCs of excluded forecast start year. Fabian.
}
if (!is.null(exp)) {
if (!obsproj) {
exp_sub <- ClimProjDiags::Subset(exp, 2, c(1:ntime)[-tt], drop = FALSE)
#NOTE: different from s2dverification. Here, 'member' is considered.
exp_EOF <- EOF(exp_sub, lat = lat, lon = lon, time_dim = names(ntime),
space_dim = c(names(nlat), names(nlon)), neofs = 1)
## Correct polarity of pattern.
#NOTE: different from s2dverification
for (imemb in 1:nmemb_exp) {
if (0 < mean(exp_EOF$EOFs[1, which.min(abs(lat - 65)), , imemb], na.rm = T)) {
exp_EOF$EOFs[1, , , imemb] <- exp_EOF$EOFs[1, , , imemb] * (-1)
}
}
# exp_EOF$PCs <- exp_EOF$PCs * sign # not used
### Lines below could be simplified further by computing
### ProjectField() only on the year of interest... (though this is
### not vital). Lauriane
PF <- ProjectField(exp, eof = exp_EOF, time_dim = names(ntime),
space_dim = c(names(nlat), names(nlon)), mode = 1)
NAOF.ver[tt, ] <- PF[tt, ]
} else {
## Project forecast anomalies on obs EOF
#NOTE: Because obs and exp have different nmemb, do ensemble mean to
# obs_EOF$EOFs first, then expand the memb dim to be the same as exp.
obs_EOF$EOFs <- apply(obs_EOF$EOFs, c(1, 2, 3), mean, na.rm = T)
obs_EOF$EOFs <- array(obs_EOF$EOFs, dim = c(dim(obs_EOF$EOFs), as.numeric(nmemb_exp)))
names(dim(obs_EOF$EOFs))[4] <- names(nmemb_obs)
PF <- ProjectField(exp, obs_EOF, mode = 1)
NAOF.ver[tt, ] <- PF[tt, ]
}
}
}
#NOTE: EOFs_obs is not returned because it's only the result of the last sdate
# (It is returned in s2dverification.)
if (!is.null(exp) & !is.null(obs)) {
return(list(exp = NAOF.ver, obs = NAOO.ver)) #, EOFs_obs = obs_EOF))
} else if (!is.null(exp)) {
return(list(exp = NAOF.ver))
} else if (!is.null(obs)) {
return(list(obs = NAOO.ver))
}
}