Newer
Older
#'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 or the forecast
#'anomalies onto the EOF pattern of the other years of the forecast.
#'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 hindcast (exp) and observations
#'(obs) based on the leading EOF pattern, or, if forecast (exp_cor) is provided,
#'the NAO index for forecast and the chosen corresponding data (exp or obs).
#'
#'@param exp A named numeric array of North Atlantic SLP (20N-80N, 80W-40E)
#' hindcast 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', '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 exp_cor A named numeric array of the Nort Atlantic SLP (20-80N, 80W-40E)
#' forecast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with
#' dimension 'time_dim' of length 1 (as in the case of an operational
#' forecast), 'memb_dim', 'ftime_dim', and 'space_dim' at least.
#' If only NAO of reference period 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', optional). If 'obs' has memb_dim, the length
#' must be 1. 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. If average is not needed, set NULL. 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. If 'exp_cor' is provided, 'obs'
#' will be used when obsproj is TRUE and 'exp' will be used when obsproj is
#' FALSE.
#'@param ncores An integer indicating the number of cores to use for parallel
#' computation. The default value is NULL.
#'
#'@return
#'A list which contains some of the following items depending on the data inputs:
#' A numeric array of hindcast NAO index in verification format with the same
#' dimensions as 'exp' except space_dim and ftime_dim. If ftime_avg is NULL,
#' ftime_dim remains.
#' A numeric array of observation NAO index in verification format with the same
#' dimensions as 'obs' except space_dim and ftime_dim. If ftime_avg is NULL,
#' ftime_dim remains.
#'\item{exp_cor}{
#' A numeric array of forecast NAO index in verification format with the same
#' dimensions as 'exp_cor' except space_dim and ftime_dim. If ftime_avg is NULL,
#' ftime_dim remains.
#' }
#'
#'@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
#'# Make up synthetic data
#'set.seed(1)
#'exp <- array(rnorm(1620), dim = c(member = 2, sdate = 3, ftime = 5, lat = 6, lon = 9))
#'set.seed(2)
#'obs <- array(rnorm(1620), dim = c(member = 1, sdate = 3, ftime = 5, lat = 6, lon = 9))
#'lat <- seq(20, 80, length.out = 6)
#'lon <- seq(-80, 40, length.out = 9)
#'nao <- NAO(exp = exp, obs = obs, lat = lat, lon = lon)
#'exp_cor <- array(rnorm(540), dim = c(member = 2, sdate = 1, ftime = 5, lat = 6, lon = 9))
#'nao <- NAO(exp = NULL, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon)
#' \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, exp_cor = 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, obs, and exp_cor (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' 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' must have at least dimensions ",
"time_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.")
}
}
if (!is.numeric(exp_cor)) {
stop("Parameter 'exp_cor' must be a numeric array.")
stop(paste0("Parameter 'exp_cor' must have at least dimensions ",
"time_dim, memb_dim, space_dim, and ftime_dim."))
}
if(any(is.null(names(dim(exp_cor)))) | any(nchar(names(dim(exp_cor))) == 0)) {
stop("Parameter 'exp_cor' 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' dimension.")
}
}
if (!is.null(obs)) {
if (!time_dim %in% names(dim(obs))) {
stop("Parameter 'time_dim' is not found in 'obs' dimension.")
if (!is.null(exp_cor)) {
if (!time_dim %in% names(dim(exp_cor))) {
stop("Parameter 'time_dim' is not found in 'exp_cor' dimension.")
}
if (dim(exp_cor)[time_dim] > 1) {
stop(paste("Parameter 'exp_cor' is expected to have length 1 on the ",
time_dim, "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' dimension.")
add_member_back <- FALSE
if (memb_dim %in% names(dim(obs))) {
if (dim(obs)[memb_dim] != 1) {
stop("The length of parameter 'memb_dim' in 'obs' must be 1.")
} else {
add_member_back <- TRUE
obs <- ClimProjDiags::Subset(obs, memb_dim, 1, drop = 'selected')
}
if (!is.null(exp_cor)) {
if (!memb_dim %in% names(dim(exp_cor))) {
stop("Parameter 'memb_dim' is not found in 'exp_cor' 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.")
}
}
if (!is.null(exp_cor)) {
if (any(!space_dim %in% names(dim(exp_cor)))) {
stop("Parameter 'space_dim' is not found in 'exp_cor' dimensions.")
}
}
## 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.")
}
}
if (!is.null(exp_cor)) {
if (!ftime_dim %in% names(dim(exp_cor))) {
stop("Parameter 'ftime_dim' is not found in 'exp_cor' dimensions.")
}
}
## 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)]
throw_error <- FALSE
if (length(name_exp) != length(name_obs)) {
throw_error <- TRUE
} else if (any(name_exp != name_obs)) {
throw_error <- TRUE
} else if(!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
throw_error <- TRUE
}
if (throw_error) {
stop(paste0("Parameter 'exp' and 'obs' must have the same names and lengths ",
"of all the dimensions except 'memb_dim'."))
}
}
## ftime_avg
if (!is.null(ftime_avg)) {
if (!is.vector(ftime_avg) | !is.numeric(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.")
}
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.")
}
if (!is.null(exp_cor)) {
if (max(ftime_avg) > dim(exp_cor)[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'."))
}
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'."))
}
}
if (!is.null(exp_cor)) {
if (!is.numeric(lat) | length(lat) != dim(exp_cor)[space_dim[1]]) {
stop(paste0("Parameter 'lat' must be a numeric vector with the same ",
"length as the latitude dimension of 'exp_cor'."))
}
if (!is.numeric(lon) | length(lon) != dim(exp_cor)[space_dim[2]]) {
stop(paste0("Parameter 'lon' must be a numeric vector with the same ",
"length as the longitude dimension of 'exp_cor'."))
}
}
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
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) & is.null(exp_cor)) {
.warning("parameter 'obsproj' set to TRUE but no 'exp' nor 'exp_cor' 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.")
}
}
if (!is.null(ftime_avg)) {
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_cor)) {
exp_cor_sub <- ClimProjDiags::Subset(exp_cor, ftime_dim, ftime_avg, drop = FALSE)
exp_cor <- MeanDims(exp_cor_sub, ftime_dim, na.rm = TRUE)
}
# wght
wght <- array(sqrt(cos(lat * pi/180)), dim = c(length(lat), length(lon)))
if (is.null(exp_cor)) {
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(time_dim, space_dim)),
fun = .NAO,
lat = lat, wght = wght,
obsproj = obsproj, add_member_back = add_member_back,
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,
lat = lat, wght = wght, obs = NULL,
obsproj = obsproj, add_member_back = FALSE,
ncores = ncores)
} else if (!is.null(obs)) {
if (add_member_back) {
output_dims <- list(obs = c(time_dim, memb_dim))
} else {
output_dims <- list(obs = time_dim)
}
res <- Apply(list(obs = obs),
target_dims = list(obs = c(time_dim, space_dim)),
output_dims = output_dims,
obsproj = obsproj, add_member_back = add_member_back,
ncores = ncores)
}
} else { # exp_cor provided
if (!is.null(exp) & !obsproj) {
.warning("Reference data provided in parameter 'obs' is not used when ",
"parameter 'obsproj' is set to FALSE and 'exp' is provided.")
res <- Apply(list(exp = exp, exp_cor = exp_cor),
target_dims = list(exp = c(memb_dim, time_dim, space_dim),
exp_cor = c(memb_dim, time_dim, space_dim)),
fun = .NAO,
obsproj = obsproj, add_member_back = add_member_back,
.warning("Experimental data provided in parameter 'exp' is not used when ",
"parameter 'obsproj' is set to TRUE and 'obs' is provided.")
res <- Apply(list(obs = obs, exp_cor = exp_cor),
target_dims = list(obs = c(time_dim, space_dim),
exp_cor = c(memb_dim, time_dim, space_dim)),
fun = .NAO,
exp = NULL,
lat = lat, wght = wght,
obsproj = obsproj, add_member_back = add_member_back,
ncores = ncores)
} else {
stop("Unconsidered case. Please review inputs 'exp', 'obs', 'exp_cor', ",
"'obsproj' or open an issue on GitLab.")
.NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, wght, obsproj = TRUE,
add_member_back = FALSE) {
# exp: [memb_exp, sdate, lat, lon]
# exp_cor: [memb, sdate = 1, lat, lon]
if (!is.null(exp)) {
ntime <- dim(exp)[2]
nlat <- dim(exp)[3]
nlon <- dim(exp)[4]
nmemb_exp <- dim(exp)[1]
} else {
ntime <- dim(obs)[1]
nlat <- dim(obs)[2]
nlon <- dim(obs)[3]
if (!is.null(exp_cor)) {
ntime_exp_cor <- dim(exp_cor)[2] # should be 1
nmemb_exp_cor <- dim(exp_cor)[1]
}
if (!is.null(obs)) NAOO.ver <- array(NA, dim = ntime)
if (!is.null(exp)) NAOF.ver <- array(NA, dim = c(ntime, nmemb_exp))
if (!is.null(exp_cor)) {
NAOF.cor <- array(NA, dim = c(ntime_exp_cor, nmemb_exp_cor))
#NOTE: The dimensions are flipped to fill in data correctly. Need to flip it back later.
if (!is.null(exp)) NAOF.ver <- array(NA, dim = c(nmemb_exp, ntime))
}
if (is.null(exp_cor)) { # cross-validation:
for (tt in 1:ntime) { #sdate
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
if (!is.null(obs)) {
## Calculate observation EOF. Excluding one forecast start year.
obs_sub <- obs[c(1:ntime)[-tt], , , drop = FALSE]
obs_EOF <- .EOF(obs_sub, neofs = 1, wght = wght) # $EOFs: [mode, lat, lon]
## Correct polarity of pattern.
# dim(obs_EOF$EOFs): [mode, lat, lon]
if (0 < mean(obs_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
obs_EOF$EOFs <- obs_EOF$EOFs * (-1)
# obs_EOF$PCs <- obs_EOF$PCs * (-1) # not used
}
## Project observed anomalies.
PF <- .ProjectField(obs, eof_mode = obs_EOF$EOFs[1, , ], wght = wght) # [sdate]
## Keep PCs of excluded forecast start year. Fabian.
NAOO.ver[tt] <- PF[tt]
}
if (!is.null(exp)) {
if (!obsproj) {
exp_sub <- exp[, c(1:ntime)[-tt], , , drop = FALSE]
# Combine 'memb' and 'sdate' to calculate EOF
dim(exp_sub) <- c(nmemb_exp * (ntime - 1), nlat, nlon)
exp_EOF <- .EOF(exp_sub, neofs = 1, wght = wght) # $EOFs: [mode, lat, lon]
## Correct polarity of pattern.
##NOTE: different from s2dverification, which doesn't use mean().
# if (0 < exp_EOF$EOFs[1, which.min(abs(lat - 65)), ]) {
if (0 < mean(exp_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
exp_EOF$EOFs <- exp_EOF$EOFs * (-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
for (imemb in 1:nmemb_exp) {
PF <- .ProjectField(exp[imemb, , , ], eof_mode = exp_EOF$EOFs[1, , ], wght = wght) # [sdate, memb]
NAOF.ver[tt, imemb] <- PF[tt]
}
} else {
## Project forecast anomalies on obs EOF
for (imemb in 1:nmemb_exp) {
PF <- .ProjectField(exp[imemb, , , ], eof_mode = obs_EOF$EOFs[1, , ], wght = wght) # [sdate]
NAOF.ver[tt, imemb] <- PF[tt]
}
}
}
} # for loop sdate
# add_member_back
if (add_member_back) {
suppressWarnings(
NAOO.ver <- InsertDim(NAOO.ver, 2, 1, name = names(dim(exp))[1])
)
}
} else { # exp_cor provided
if (!is.null(obs) & obsproj) {
obs_EOF <- .EOF(obs, neofs = 1, wght = wght) # $EOFs: [mode, lat, lon]
## Correct polarity of pattern.
# dim(obs_EOF$EOFs): [mode, lat, lon]
if (0 < mean(obs_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
obs_EOF$EOFs <- obs_EOF$EOFs * (-1)
PF <- .ProjectField(obs, eof_mode = obs_EOF$EOFs[1, , ], wght = wght) # [sdate]
NAOO.ver[] <- PF
## Project observed anomalies.
for (imemb in 1:dim(exp_cor)[1]) {
exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb,
PF <- .ProjectField(exp_sub, eof_mode = obs_EOF$EOFs[1, , ], wght = wght) # [sdate]
NAOF.cor[imemb] <- PF
} else if (!is.null(exp) & !obsproj) {
dim(exp) <- c(nmemb_exp*ntime, nlat, nlon)
exp_EOF <- .EOF(exp, neofs = 1, wght = wght) # $EOFs: [mode, lat, lon]
## Correct polarity of pattern.
# dim(obs_EOF$EOFs): [mode, lat, lon]
if (0 < mean(exp_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
exp_EOF$EOFs <- exp_EOF$EOFs * (-1)
# obs_EOF$PCs <- obs_EOF$PCs * (-1) # not used
}
PF <- .ProjectField(exp, eof_mode = exp_EOF$EOFs[1, , ], wght = wght) # [sdate]
NAOF.ver[, ] <- PF
# Flip the dimensions back
NAOF.ver <- aperm(NAOF.ver, 2:1)
## Project observed anomalies.
for (imemb in 1:dim(exp_cor)[1]) {
exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb,
PF <- .ProjectField(exp_sub, eof_mode = exp_EOF$EOFs[1, , ], wght = wght) # [sdate]
NAOF.cor[imemb] <- PF
}
}
#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_cor)) {
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))
}
} else {
if (!is.null(exp) & !obsproj) {
return(list(exp = NAOF.ver, exp_cor = NAOF.cor))
} else if (!is.null(obs) & obsproj) {
return(list(obs = NAOO.ver, exp_cor = NAOF.cor))
}