From 6c3bda23ef59db869e2c1404d39c194972a799fb Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 30 Aug 2023 11:45:29 +0200 Subject: [PATCH 01/12] First attemp forecast NAO --- R/NAO.R | 245 +++++++++++++++++++++++++++++++++++------------ R/ProjectField.R | 1 - 2 files changed, 182 insertions(+), 64 deletions(-) diff --git a/R/NAO.R b/R/NAO.R index fb5220c..60286f4 100644 --- a/R/NAO.R +++ b/R/NAO.R @@ -11,7 +11,7 @@ #'(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 +#' 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. @@ -20,6 +20,10 @@ #' 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 \cod{Ano()} or \code{Ano_CrossValid()] with +#' If only NAO of observational 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 @@ -73,6 +77,8 @@ #'lon <- seq(-80, 40, length.out = 9) #'nao <- NAO(exp = exp, obs = obs, lat = lat, lon = lon) #' +#'exp_cor <- array(rnorm(2*5*6*9), dim = c(member = 2, sdate = 1, ftime = 5, lat = 6, lon = 9)) +#'nao <- NAO(exp = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon) #'# plot the NAO index #' \dontrun{ #'nao$exp <- Reorder(nao$exp, c(2, 1)) @@ -84,7 +90,7 @@ #'@import multiApply #'@importFrom ClimProjDiags Subset #'@export -NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate', +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) { @@ -118,6 +124,18 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate', stop("Parameter 'obs' must have dimension names.") } } + if (!is.null(exp_cor)) { + if (!is.numeric(exp)) { + stop("Parameter 'exp' must be a numeric array.") + } + if (is.null(dim(exp_cor))) { + 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_cor)))) | any(nchar(names(dim(exp_cor))) == 0)) { + stop("Parameter 'exp' must have dimension names.") + } + } ## time_dim if (!is.character(time_dim) | length(time_dim) > 1) { stop("Parameter 'time_dim' must be a character string.") @@ -153,6 +171,11 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate', add_member_back <- FALSE } } + 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.") @@ -167,6 +190,11 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate', 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.") @@ -181,6 +209,11 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate', 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))) @@ -208,11 +241,17 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate', 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 (!is.null(obs)) { 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)) { @@ -234,7 +273,8 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate', stop(paste0("Parameter 'lon' must be a numeric vector with the same ", "length as the longitude dimension of 'exp' and 'obs'.")) } - } else { + } + if (!is.null(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'.")) @@ -244,6 +284,16 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate', "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'.")) + } + } stop_needed <- FALSE if (max(lat) > 80 | min(lat) < 20) { stop_needed <- TRUE @@ -303,19 +353,34 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate', 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) & !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) + if (is.null(exp_cor)) { + 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 { + res <- Apply(list(exp, obs, exp_cor), + target_dims = list(exp = c(memb_dim, time_dim, space_dim), + obs = c(time_dim, space_dim), + exp_cor = c(memb_dim, 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)), @@ -340,11 +405,12 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate', return(res) } -.NAO <- function(exp = NULL, obs = NULL, lat, wght, obsproj = TRUE, add_member_back = FALSE) { +.NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, + lat, wght, obsproj = TRUE, add_member_back = FALSE) { # exp: [memb_exp, sdate, lat, lon] # obs: [sdate, lat, lon] + # exp_cor: [memb, sdate=1, lat, lon] # wght: [lat, lon] - if (!is.null(exp)) { ntime <- dim(exp)[2] nlat <- dim(exp)[3] @@ -358,73 +424,126 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate', 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 = dim(exp_cor)[1]) - for (tt in 1:ntime) { #sdate + if (is.null(exp_cor)) { # cross-validation: + for (tt in 1:ntime) { #sdate - 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] + 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) -# obs_EOF$PCs <- obs_EOF$PCs * (-1) # not used +# 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] + NAOO.ver <- PF - ## 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 - } + ## Project observed anomalies. + for (imemb in 1:dim(exp_cor)[1]) { + exp_sub <- ClimProjDiags::Subset(exp_cor, along = memb_dim, indices = imemb, + drop = 'selected') - ### 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] - } + 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 - } # for loop sdate + ## Project observed anomalies. + for (imemb in 1:dim(exp_cor)[1]) { + exp_sub <- ClimProjDiags::Subset(exp_cor, along = memb_dim, indices = imemb, + drop = 'selected') - # add_member_back - if (add_member_back) { - suppressWarnings( - NAOO.ver <- InsertDim(NAOO.ver, 2, 1, name = names(dim(exp))[1]) - ) + 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) & !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)) + 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)) + } else { + stop("Unconsidered case. Please review") + } } } diff --git a/R/ProjectField.R b/R/ProjectField.R index 55e7fd2..7e73253 100644 --- a/R/ProjectField.R +++ b/R/ProjectField.R @@ -229,7 +229,6 @@ ProjectField <- function(ano, eof, time_dim = 'sdate', space_dim = c('lat', 'lon # wght: [lat, lon] ntime <- dim(ano)[1] - if (length(dim(eof_mode)) == 2) { # mode != NULL # Initialization of pc.ver. pc.ver <- array(NA, dim = ntime) #[sdate] -- GitLab From ccac96527b80568ac9faf219ccbf55476fadc512 Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 30 Aug 2023 15:53:01 +0200 Subject: [PATCH 02/12] exp_cor checks correctly placed --- R/NAO.R | 83 ++++++++++++++++++++++++++++++++------------------------- 1 file changed, 46 insertions(+), 37 deletions(-) diff --git a/R/NAO.R b/R/NAO.R index 60286f4..24f8ae3 100644 --- a/R/NAO.R +++ b/R/NAO.R @@ -125,15 +125,15 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd } } if (!is.null(exp_cor)) { - if (!is.numeric(exp)) { - stop("Parameter 'exp' must be a numeric array.") + if (!is.numeric(exp_cor)) { + stop("Parameter 'exp_cor' must be a numeric array.") } if (is.null(dim(exp_cor))) { - stop(paste0("Parameter 'exp' must have at least dimensions ", + 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' must have dimension names.") + stop("Parameter 'exp_cor' must have dimension names.") } } ## time_dim @@ -329,8 +329,8 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd 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.") + if (is.null(exp) & is.null(exp_cor)) { + .warning("parameter 'obsproj' set to TRUE but no 'exp' nor 'exp_cor' provided.") } } ## ncores @@ -361,46 +361,56 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd # wght wght <- array(sqrt(cos(lat * pi/180)), dim = c(length(lat), length(lon))) - - if (!is.null(exp) & !is.null(obs)) { - if (is.null(exp_cor)) { - res <- Apply(list(exp, obs), - target_dims = list(exp = c(memb_dim, time_dim, space_dim), - obs = c(time_dim, space_dim)), + 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, fun = .NAO, - lat = lat, wght = wght, + lat = lat, wght = wght, exp = NULL, obsproj = obsproj, add_member_back = add_member_back, ncores = ncores) - } else { - res <- Apply(list(exp, obs, exp_cor), + } + } else { # exp_cor provided + if (!is.null(exp) & !obsproj) { + res <- Apply(list(exp = exp, exp_cor = exp_cor), target_dims = list(exp = c(memb_dim, time_dim, space_dim), - obs = c(time_dim, space_dim), exp_cor = c(memb_dim, time_dim, space_dim)), fun = .NAO, + obs = NULL, + lat = lat, wght = wght, + obsproj = obsproj, add_member_back = add_member_back, + ncores = ncores) + } else if (!is.null(obs) & obsproj) { + 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 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, - fun = .NAO, - lat = lat, wght = wght, exp = NULL, - obsproj = obsproj, add_member_back = add_member_back, - ncores = ncores) } return(res) } @@ -515,7 +525,6 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd } PF <- .ProjectField(exp, eof_mode = exp_EOF$EOFs[1, , ], wght = wght) # [sdate] NAOF.ver <- PF - ## Project observed anomalies. for (imemb in 1:dim(exp_cor)[1]) { exp_sub <- ClimProjDiags::Subset(exp_cor, along = memb_dim, indices = imemb, -- GitLab From e31a37a7928d870e92fb1e44e5a72c7a26434a87 Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 30 Aug 2023 16:26:34 +0200 Subject: [PATCH 03/12] Fix for only exp provided --- R/NAO.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/NAO.R b/R/NAO.R index 24f8ae3..16fd4e8 100644 --- a/R/NAO.R +++ b/R/NAO.R @@ -399,7 +399,7 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd fun = .NAO, obs = NULL, lat = lat, wght = wght, - obsproj = obsproj, add_member_back = add_member_back, + obsproj = obsproj, #add_member_back = add_member_back, ncores = ncores) } else if (!is.null(obs) & obsproj) { res <- Apply(list(obs = obs, exp_cor = exp_cor), -- GitLab From 66e6a2bae45f8433bf639d133619791c7a460f80 Mon Sep 17 00:00:00 2001 From: aho Date: Fri, 8 Sep 2023 13:05:59 +0200 Subject: [PATCH 04/12] Trivial adjustment and fix unit test --- R/NAO.R | 6 +++--- man/NAO.Rd | 10 +++++++++- tests/testthat/test-NAO.R | 2 +- 3 files changed, 13 insertions(+), 5 deletions(-) diff --git a/R/NAO.R b/R/NAO.R index 16fd4e8..a34584a 100644 --- a/R/NAO.R +++ b/R/NAO.R @@ -21,7 +21,7 @@ #' 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 \cod{Ano()} or \code{Ano_CrossValid()] with +#' forecast anomalies from \cod{Ano()} or \code{Ano_CrossValid()} with #' If only NAO of observational 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'. @@ -533,6 +533,8 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd PF <- .ProjectField(exp_sub, eof_mode = exp_EOF$EOFs[1, , ], wght = wght) # [sdate] NAOF.cor[imemb] <- PF } + } else { + stop("Unconsidered case. Please review.") } } @@ -551,8 +553,6 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd return(list(exp = NAOF.ver, exp_cor = NAOF.cor)) } else if (!is.null(obs) & obsproj) { return(list(obs = NAOO.ver, exp_cor = NAOF.cor)) - } else { - stop("Unconsidered case. Please review") } } } diff --git a/man/NAO.Rd b/man/NAO.Rd index 999fd75..e605469 100644 --- a/man/NAO.Rd +++ b/man/NAO.Rd @@ -7,6 +7,7 @@ NAO( exp = NULL, obs = NULL, + exp_cor = NULL, lat, lon, time_dim = "sdate", @@ -20,7 +21,7 @@ NAO( } \arguments{ \item{exp}{A named numeric array of North Atlantic SLP (20N-80N, 80W-40E) -forecast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with +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.} @@ -31,6 +32,11 @@ 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.} +\item{exp_cor}{A named numeric array of the Nort Atlantic SLP (20-80N, 80W-40E) +forecast anomalies from \cod{Ano()} or \code{Ano_CrossValid()} with +If only NAO of observational data needs to be computed, this parameter can +be left to NULL. The default value is NULL.} + \item{lat}{A vector of the latitudes of 'exp' and 'obs'.} \item{lon}{A vector of the longitudes of 'exp' and 'obs'.} @@ -97,6 +103,8 @@ 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(2*5*6*9), dim = c(member = 2, sdate = 1, ftime = 5, lat = 6, lon = 9)) +nao <- NAO(exp = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon) # plot the NAO index \dontrun{ nao$exp <- Reorder(nao$exp, c(2, 1)) diff --git a/tests/testthat/test-NAO.R b/tests/testthat/test-NAO.R index 91b9943..bbc7832 100644 --- a/tests/testthat/test-NAO.R +++ b/tests/testthat/test-NAO.R @@ -153,7 +153,7 @@ test_that("1. Input checks", { ) # ncores expect_error( - NAO(exp1, obs1, lat1, lon1, ncore = 3.5), + NAO(exp1, obs1, lat = lat1, lon = lon1, ncore = 3.5), "Parameter 'ncores' must be a positive integer." ) -- GitLab From 59260099aa72c4ac5d6d25314281043d99c8b36f Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 20 Sep 2023 11:57:29 +0200 Subject: [PATCH 05/12] Add missing sentence in doc --- R/NAO.R | 3 ++- tests/testthat/test-NAO.R | 8 +++++++- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/R/NAO.R b/R/NAO.R index a34584a..0d9f9a5 100644 --- a/R/NAO.R +++ b/R/NAO.R @@ -21,7 +21,8 @@ #' 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 \cod{Ano()} or \code{Ano_CrossValid()} with +#' forecast anomalies from \cod{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 lat A vector of the latitudes of 'exp' and 'obs'. diff --git a/tests/testthat/test-NAO.R b/tests/testthat/test-NAO.R index bbc7832..f8e2b1f 100644 --- a/tests/testthat/test-NAO.R +++ b/tests/testthat/test-NAO.R @@ -14,7 +14,8 @@ obs2 <- array(rnorm(108), dim = c(sdate = 3, ftime = 4, lat = 3, lon = 3)) lat2 <- c(80, 50, 20) lon2 <- c(-80, 0, 40) - + set.seed(3) + exp2_cor <- array(rnorm(72), dim = c(sdate = 1, ftime = 4, member = 2, lat = 3, lon = 3)) ############################################## test_that("1. Input checks", { @@ -253,6 +254,11 @@ test_that("3. dat2", { tolerance = 0.00001 ) + # exp_cor +# expect_equal( +# dim(NAO(exp2, obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp), +# c(sdate = 3, member = 2) +# ) }) ############################################## -- GitLab From 9db4e45b4b48eb6f0833f0fff3521b5157e957d0 Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 20 Sep 2023 17:39:18 +0200 Subject: [PATCH 06/12] update docu and checks --- R/NAO.R | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/R/NAO.R b/R/NAO.R index a34584a..de2454e 100644 --- a/R/NAO.R +++ b/R/NAO.R @@ -22,8 +22,9 @@ #' 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 \cod{Ano()} or \code{Ano_CrossValid()} with -#' If only NAO of observational data needs to be computed, this parameter can -#' be left to NULL. The default value is NULL. +#' dimension 'time_dim' of length 1 (as in the case of an operational forecast). +#' 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 @@ -44,7 +45,8 @@ #' 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. +#' 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. #' @@ -135,6 +137,10 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd if(any(is.null(names(dim(exp_cor)))) | any(nchar(names(dim(exp_cor))) == 0)) { stop("Parameter 'exp_cor' must have dimension names.") } + if (dim(exp_cor)[time_dim] > 1) { + stop(paste("Parameter 'exp_cor' is expected to have length 1 on the", + time_dim, "dimension.")) + } } ## time_dim if (!is.character(time_dim) | length(time_dim) > 1) { @@ -393,6 +399,10 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd } } else { # exp_cor provided if (!is.null(exp) & !obsproj) { + if (!is.null(obs)) { + 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)), @@ -402,6 +412,10 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd obsproj = obsproj, #add_member_back = add_member_back, ncores = ncores) } else if (!is.null(obs) & obsproj) { + if (!is.null(exp)) { + 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)), -- GitLab From 972d06e5f1d821450b637102e8a610086f2eab15 Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 21 Sep 2023 14:47:29 +0200 Subject: [PATCH 07/12] fix memb_dim --- R/NAO.R | 6 +++--- man/NAO.Rd | 8 +++++--- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/R/NAO.R b/R/NAO.R index de2454e..9dd0854 100644 --- a/R/NAO.R +++ b/R/NAO.R @@ -21,7 +21,7 @@ #' 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 \cod{Ano()} or \code{Ano_CrossValid()} with +#' forecast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with #' dimension 'time_dim' of length 1 (as in the case of an operational forecast). #' If only NAO of reference period needs to be computed, this parameter can #' be left to NULL. The default value is NULL. @@ -522,7 +522,7 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd ## Project observed anomalies. for (imemb in 1:dim(exp_cor)[1]) { - exp_sub <- ClimProjDiags::Subset(exp_cor, along = memb_dim, indices = imemb, + exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb, drop = 'selected') PF <- .ProjectField(exp_sub, eof_mode = obs_EOF$EOFs[1, , ], wght = wght) # [sdate] @@ -541,7 +541,7 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd NAOF.ver <- PF ## Project observed anomalies. for (imemb in 1:dim(exp_cor)[1]) { - exp_sub <- ClimProjDiags::Subset(exp_cor, along = memb_dim, indices = imemb, + exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb, drop = 'selected') PF <- .ProjectField(exp_sub, eof_mode = exp_EOF$EOFs[1, , ], wght = wght) # [sdate] diff --git a/man/NAO.Rd b/man/NAO.Rd index e605469..84f445e 100644 --- a/man/NAO.Rd +++ b/man/NAO.Rd @@ -33,8 +33,9 @@ If only NAO of experimental data needs to be computed, this parameter can be left to NULL. The default value is NULL.} \item{exp_cor}{A named numeric array of the Nort Atlantic SLP (20-80N, 80W-40E) -forecast anomalies from \cod{Ano()} or \code{Ano_CrossValid()} with -If only NAO of observational data needs to be computed, this parameter can +forecast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with +dimension 'time_dim' of length 1 (as in the case of an operational forecast). +If only NAO of reference period needs to be computed, this parameter can be left to NULL. The default value is NULL.} \item{lat}{A vector of the latitudes of 'exp' and 'obs'.} @@ -64,7 +65,8 @@ 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.} +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.} \item{ncores}{An integer indicating the number of cores to use for parallel computation. The default value is NULL.} -- GitLab From 874f6d25d04f54bef0c74eaad42c8d9321ce5c93 Mon Sep 17 00:00:00 2001 From: aho Date: Fri, 22 Sep 2023 16:20:33 +0200 Subject: [PATCH 08/12] Correct exp output dim when exp and exp_cor are used --- R/NAO.R | 12 +++++++++--- tests/testthat/test-NAO.R | 32 ++++++++++++++++++++++++++++---- 2 files changed, 37 insertions(+), 7 deletions(-) diff --git a/R/NAO.R b/R/NAO.R index 40fc533..0f3cb7f 100644 --- a/R/NAO.R +++ b/R/NAO.R @@ -470,8 +470,11 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd 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)) - + 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 @@ -561,7 +564,10 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd # obs_EOF$PCs <- obs_EOF$PCs * (-1) # not used } PF <- .ProjectField(exp, eof_mode = exp_EOF$EOFs[1, , ], wght = wght) # [sdate] - NAOF.ver <- PF + 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, diff --git a/tests/testthat/test-NAO.R b/tests/testthat/test-NAO.R index 5f1eafa..7b64cac 100644 --- a/tests/testthat/test-NAO.R +++ b/tests/testthat/test-NAO.R @@ -263,10 +263,34 @@ test_that("3. dat2", { dim(NAO(exp = NULL, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$obs), c(sdate = 3) ) -#res <- NAO(exp = NULL, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2, ftime_avg = NULL) -#WRONG -#res1 <- NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F) -#res2 <- NAO(exp = exp2, obs = obs2, exp_cor = NULL, lat = lat2, lon = lon2) + expect_equal( + c((NAO(exp = NULL, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp_cor)), + c(0.2121340, 0.1634516), + tolerance = 0.0001) + expect_equal( + c((NAO(exp = NULL, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$obs)), + c(0.3511294, -0.7196260, -1.5123894), + tolerance = 0.0001 + ) + + expect_equal( + dim(NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)$exp_cor), + c(sdate = 1, member = 2) + ) + expect_equal( + dim(NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)$exp), + c(sdate = 3, member = 2) + ) + expect_equal( + c(NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)$exp), + c(-0.4225802, -1.3616918, 0.2808729, 0.2723077, -0.7584804, -0.7264514), + tolerance = 0.0001 + ) + expect_equal( + c(NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)$exp_cor), + c(-0.22949531, -0.06946422), + tolerance = 0.0001 + ) }) -- GitLab From d00a030a9db883a9fb65cbc6c994bfe43111b4d6 Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 24 Oct 2023 17:40:43 +0200 Subject: [PATCH 09/12] NAO fcst and hcst returned --- R/NAO.R | 88 +++++++++++++++++++++++---------------------------------- 1 file changed, 36 insertions(+), 52 deletions(-) diff --git a/R/NAO.R b/R/NAO.R index 0f3cb7f..33adc47 100644 --- a/R/NAO.R +++ b/R/NAO.R @@ -88,7 +88,7 @@ #'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) +#'nao <- NAO(exp = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon, obsproj = T) #'# plot the NAO index #' \dontrun{ #'nao$exp <- Reorder(nao$exp, c(2, 1)) @@ -104,7 +104,6 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd 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)) { @@ -411,36 +410,17 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd ncores = ncores) } } else { # exp_cor provided - if (!is.null(exp) & !obsproj) { - if (!is.null(obs)) { - .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), + if (!is.null(exp) & !is.null(obs) & !is.null(exp_cor)) { + res <- Apply(list(exp = exp, obs = obs, exp_cor = exp_cor), target_dims = list(exp = c(memb_dim, time_dim, space_dim), + obs = c(time_dim, space_dim), exp_cor = c(memb_dim, time_dim, space_dim)), fun = .NAO, - obs = NULL, lat = lat, wght = wght, obsproj = obsproj, add_member_back = add_member_back, ncores = ncores) - - } else if (!is.null(obs) & obsproj) { - if (!is.null(exp)) { - .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.") + } else { + stop("Parameters 'exp' and 'obs' are required when 'exp_cor' is not null.") } } @@ -467,13 +447,12 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd 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)) NAOF.ver <- array(NA, dim = c(nmemb_exp, ntime)) } if (is.null(exp_cor)) { # cross-validation: for (tt in 1:ntime) { #sdate @@ -535,26 +514,32 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd } } 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) -# obs_EOF$PCs <- obs_EOF$PCs * (-1) # not used - } - PF <- .ProjectField(obs, eof_mode = obs_EOF$EOFs[1, , ], wght = wght) # [sdate] - NAOO.ver[] <- PF - - ## Project observed anomalies. + ## Calculate observation EOF. Without cross-validation + 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) + # obs_EOF$PCs <- obs_EOF$PCs * (-1) # not used + } + ## Project observed anomalies. + PF <- .ProjectField(obs, eof_mode = obs_EOF$EOFs, wght = wght) # [sdate] + NAOO.ver[] <- PF[1,] + if (obsproj) { + ## Project forecasts anomalies for (imemb in 1:dim(exp_cor)[1]) { exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb, drop = 'selected') - PF <- .ProjectField(exp_sub, eof_mode = obs_EOF$EOFs[1, , ], wght = wght) # [sdate] - NAOF.cor[imemb] <- PF + NAOF.cor[,imemb] <- PF } - } else if (!is.null(exp) & !obsproj) { + for (imemb in 1:dim(exp)[1]) { + exp_sub <- ClimProjDiags::Subset(exp, along = 1, indices = imemb, + drop = 'selected') + PF <- .ProjectField(exp_sub, eof_mode = obs_EOF$EOFs[1, , ], wght = wght) + NAOF.ver[,imemb] <- PF + } + } else if (!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. @@ -563,12 +548,16 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd 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 + dim(exp) <- c(nmemb_exp, ntime, nlat, nlon) + for (imemb in 1:dim(exp)[1]) { + exp_sub <- ClimProjDiags::Subset(exp, along = 1, indices = imemb, drop = 'selected') + PF <- .ProjectField(exp_sub, eof_mode = exp_EOF$EOFs[1, , ], wght = wght) # [sdate] + NAOF.ver[ ,imemb] <- PF + } # Flip the dimensions back NAOF.ver <- aperm(NAOF.ver, 2:1) - ## Project observed anomalies. + ## Project hindcast anomalies on forecast for (imemb in 1:dim(exp_cor)[1]) { exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb, drop = 'selected') @@ -578,7 +567,6 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd } } } - #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)) { @@ -590,10 +578,6 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd 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)) - } + return(list(exp = NAOF.ver, obs = NAOO.ver, exp_cor = NAOF.cor)) } } -- GitLab From 600907e9aae46875cf6eb37d7432203e357be2a5 Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 15 Nov 2023 14:31:29 +0100 Subject: [PATCH 10/12] Tidy codes, fix unit tests --- R/NAO.R | 214 ++++++++++++++++++-------------------- man/NAO.Rd | 16 +-- tests/testthat/test-NAO.R | 73 +++++++++---- 3 files changed, 161 insertions(+), 142 deletions(-) diff --git a/R/NAO.R b/R/NAO.R index 33adc47..2a0e118 100644 --- a/R/NAO.R +++ b/R/NAO.R @@ -9,7 +9,7 @@ #'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). +#'the NAO index for forecast and the corresponding data (exp and 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 @@ -44,12 +44,12 @@ #' 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. +#' reference (TRUE, default) or compute the NAO by first computing the leading +#' EOF of the forecast anomalies (in cross-validation mode, i.e. leave the +#' evaluated year out), then projecting forecast anomalies onto this EOF +#' (FALSE). If 'exp_cor' is provided, 'obs' will be used when obsproj is TRUE +#' and 'exp' will be used when obsproj is FALSE, and no cross-validation is +#' applied. #'@param ncores An integer indicating the number of cores to use for parallel #' computation. The default value is NULL. #' @@ -141,9 +141,12 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd 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)) { + if (any(is.null(names(dim(exp_cor)))) | any(nchar(names(dim(exp_cor))) == 0)) { stop("Parameter 'exp_cor' must have dimension names.") } + if (is.null(exp) || is.null(obs)) { + stop("Parameters 'exp' and 'obs' are required when 'exp_cor' is not provided.") + } } ## time_dim if (!is.character(time_dim) | length(time_dim) > 1) { @@ -164,8 +167,8 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd 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.")) + stop("Parameter 'exp_cor' is expected to have length 1 in ", + time_dim, "dimension.") } } @@ -233,6 +236,7 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd } } ## exp and obs (2) + #TODO: Add checks for exp_cor if (!is.null(exp) & !is.null(obs)) { name_exp <- sort(names(dim(exp))) name_obs <- sort(names(dim(obs))) @@ -242,12 +246,12 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd throw_error <- TRUE } else if (any(name_exp != name_obs)) { throw_error <- TRUE - } else if(!all(dim(exp)[name_exp] == dim(obs)[name_obs])) { + } 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'.")) + stop("Parameter 'exp' and 'obs' must have the same names and lengths ", + "of all the dimensions except 'memb_dim'.") } } ## ftime_avg @@ -410,18 +414,14 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd ncores = ncores) } } else { # exp_cor provided - if (!is.null(exp) & !is.null(obs) & !is.null(exp_cor)) { - res <- Apply(list(exp = exp, obs = obs, exp_cor = exp_cor), - target_dims = list(exp = c(memb_dim, time_dim, space_dim), - obs = c(time_dim, space_dim), - exp_cor = c(memb_dim, time_dim, space_dim)), - fun = .NAO, - lat = lat, wght = wght, - obsproj = obsproj, add_member_back = add_member_back, - ncores = ncores) - } else { - stop("Parameters 'exp' and 'obs' are required when 'exp_cor' is not null.") - } + res <- Apply(list(exp = exp, obs = obs, exp_cor = exp_cor), + target_dims = list(exp = c(memb_dim, time_dim, space_dim), + obs = c(time_dim, space_dim), + exp_cor = c(memb_dim, time_dim, space_dim)), + fun = .NAO, + lat = lat, wght = wght, + obsproj = obsproj, add_member_back = add_member_back, + ncores = ncores) } return(res) @@ -433,6 +433,7 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd # obs: [sdate, lat, lon] # exp_cor: [memb, sdate = 1, lat, lon] # wght: [lat, lon] + if (!is.null(exp)) { ntime <- dim(exp)[2] nlat <- dim(exp)[3] @@ -447,30 +448,31 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd 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(obs)) nao_obs <- array(NA, dim = ntime) + if (!is.null(exp)) nao_exp <- array(NA, dim = c(ntime, nmemb_exp)) if (!is.null(exp_cor)) { - NAOF.cor <- array(NA, dim = c(ntime_exp_cor, nmemb_exp_cor)) + nao_exp_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 + + if (is.null(exp_cor)) { + + for (tt in 1:ntime) { # cross-validation 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 + EOF_obs <- .EOF(obs_sub, neofs = 1, wght = wght)$EOFs # [mode = 1, lat, lon] + ## Correct polarity of pattern + # EOF_obs: [mode = 1, lat, lon] + if (0 < mean(EOF_obs[1, which.min(abs(lat - 65)), ], na.rm = T)) { + EOF_obs <- EOF_obs * (-1) } ## Project observed anomalies. - PF <- .ProjectField(obs, eof_mode = obs_EOF$EOFs[1, , ], wght = wght) # [sdate] + PF <- .ProjectField(obs, eof_mode = EOF_obs[1, , ], wght = wght) # [sdate] ## Keep PCs of excluded forecast start year. Fabian. - NAOO.ver[tt] <- PF[tt] + nao_obs[tt] <- PF[tt] } if (!is.null(exp)) { @@ -478,106 +480,94 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd 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] + EOF_exp <- .EOF(exp_sub, neofs = 1, wght = wght)$EOFs # [mode = 1, lat, lon] - ## Correct polarity of pattern. + ## 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 +# if (0 < EOF_exp[1, which.min(abs(lat - 65)), ]) { + if (0 < mean(EOF_exp[1, which.min(abs(lat - 65)), ], na.rm = T)) { + EOF_exp <- EOF_exp * (-1) } ### 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] + PF <- .ProjectField(exp[imemb, , , ], eof_mode = EOF_exp[1, , ], wght = wght) # [sdate, memb] + nao_exp[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] + PF <- .ProjectField(exp[imemb, , , ], eof_mode = EOF_obs[1, , ], wght = wght) # [sdate] + nao_exp[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 + ## Calculate observation EOF. Without cross-validation - 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) - # obs_EOF$PCs <- obs_EOF$PCs * (-1) # not used - } - ## Project observed anomalies. - PF <- .ProjectField(obs, eof_mode = obs_EOF$EOFs, wght = wght) # [sdate] - NAOO.ver[] <- PF[1,] - if (obsproj) { - ## Project forecasts anomalies - for (imemb in 1:dim(exp_cor)[1]) { - exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb, - drop = 'selected') - PF <- .ProjectField(exp_sub, eof_mode = obs_EOF$EOFs[1, , ], wght = wght) # [sdate] - NAOF.cor[,imemb] <- PF - } - for (imemb in 1:dim(exp)[1]) { - exp_sub <- ClimProjDiags::Subset(exp, along = 1, indices = imemb, - drop = 'selected') - PF <- .ProjectField(exp_sub, eof_mode = obs_EOF$EOFs[1, , ], wght = wght) - NAOF.ver[,imemb] <- PF - } - } else if (!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 - } - dim(exp) <- c(nmemb_exp, ntime, nlat, nlon) - for (imemb in 1:dim(exp)[1]) { - exp_sub <- ClimProjDiags::Subset(exp, along = 1, indices = imemb, drop = 'selected') - PF <- .ProjectField(exp_sub, eof_mode = exp_EOF$EOFs[1, , ], wght = wght) # [sdate] - NAOF.ver[ ,imemb] <- PF + EOF_obs <- .EOF(obs, neofs = 1, wght = wght)$EOFs # [mode = 1, lat, lon] + ## Correct polarity of pattern + # EOF_obs: [mode, lat, lon] + if (0 < mean(EOF_obs[1, which.min(abs(lat - 65)), ], na.rm = T)) { + EOF_obs <- EOF_obs * (-1) + } + ## Project observed anomalies + PF <- .ProjectField(obs, eof_mode = EOF_obs, wght = wght) # [mode = 1, sdate] + nao_obs[] <- PF[1, ] + + if (!obsproj) { + # Calculate EOF_exp + tmp <- array(exp, dim = c(nmemb_exp * ntime, nlat, nlon)) + EOF_exp <- .EOF(tmp, neofs = 1, wght = wght)$EOFs # [mode = 1, lat, lon] + ## Correct polarity of pattern + if (0 < mean(EOF_exp[1, which.min(abs(lat - 65)), ], na.rm = T)) { + EOF_exp <- EOF_exp * (-1) } - # Flip the dimensions back - NAOF.ver <- aperm(NAOF.ver, 2:1) + eof_mode_input <- EOF_exp[1, , ] + } else { + eof_mode_input <- EOF_obs[1, , ] + } - ## Project hindcast anomalies on forecast - for (imemb in 1:dim(exp_cor)[1]) { - exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb, - drop = 'selected') + # Calculate NAO_exp + for (imemb in 1:dim(exp)[1]) { + exp_sub <- ClimProjDiags::Subset(exp, along = 1, indices = imemb, + drop = 'selected') + PF <- .ProjectField(exp_sub, eof_mode = eof_mode_input, wght = wght) # [sdate] + nao_exp[ , imemb] <- PF + } - PF <- .ProjectField(exp_sub, eof_mode = exp_EOF$EOFs[1, , ], wght = wght) # [sdate] - NAOF.cor[imemb] <- PF - } + # Calculate NAO_exp_cor + for (imemb in 1:dim(exp_cor)[1]) { + exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb, + drop = 'selected') + PF <- .ProjectField(exp_sub, eof_mode = eof_mode_input, wght = wght) # [sdate] + nao_exp_cor[, imemb] <- PF } + } - #NOTE: EOFs_obs is not returned because it's only the result of the last sdate - # (It is returned in s2dverification.) + # add_member_back + if (add_member_back) { + memb_dim_name <- ifelse(!is.null(names(dim(exp))[1]), names(dim(exp))[1], 'member') + nao_obs <- InsertDim(nao_obs, 2, 1, name = memb_dim_name) + } + + # Return results 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)) + res <- NULL + if (!is.null(exp)) { + res <- c(res, list(exp = nao_exp)) } + if (!is.null(obs)) { + res <- c(res, list(obs = nao_obs)) + } + return(res) + } else { - return(list(exp = NAOF.ver, obs = NAOO.ver, exp_cor = NAOF.cor)) + return(list(exp = nao_exp, obs = nao_obs, exp_cor = nao_exp_cor)) } } diff --git a/man/NAO.Rd b/man/NAO.Rd index c8ab94a..8545e6b 100644 --- a/man/NAO.Rd +++ b/man/NAO.Rd @@ -63,12 +63,12 @@ value is 2:4, i.e., from 2nd to 4th forecast time steps.} \item{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.} +reference (TRUE, default) or compute the NAO by first computing the leading +EOF of the forecast anomalies (in cross-validation mode, i.e. leave the +evaluated year out), then projecting forecast anomalies onto this EOF +(FALSE). If 'exp_cor' is provided, 'obs' will be used when obsproj is TRUE +and 'exp' will be used when obsproj is FALSE, and no cross-validation is +applied.} \item{ncores}{An integer indicating the number of cores to use for parallel computation. The default value is NULL.} @@ -101,7 +101,7 @@ 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). +the NAO index for forecast and the corresponding data (exp and obs). } \examples{ # Make up synthetic data @@ -114,7 +114,7 @@ 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) +nao <- NAO(exp = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon, obsproj = T) # plot the NAO index \dontrun{ nao$exp <- Reorder(nao$exp, c(2, 1)) diff --git a/tests/testthat/test-NAO.R b/tests/testthat/test-NAO.R index 7b64cac..c3f0f75 100644 --- a/tests/testthat/test-NAO.R +++ b/tests/testthat/test-NAO.R @@ -6,7 +6,9 @@ obs1 <- array(rnorm(72), dim = c(dataset = 1, member = 1, sdate = 3, ftime = 4, lat = 2, lon = 3)) lat1 <- c(20, 80) lon1 <- c(40, 280, 350) - + set.seed(3) + exp1_cor <- array(rnorm(72), dim = c(sdate = 1, ftime = 4, member = 3, lat = 2, lon = 3)) + # dat2 set.seed(1) exp2 <- array(rnorm(216), dim = c(sdate = 3, ftime = 4, member = 2, lat = 3, lon = 3)) @@ -221,6 +223,38 @@ test_that("2. dat1", { c(NAO(exp1, obs1, lat = lat1, lon = lon1, ftime_avg = 1)$exp) ) + # exp_cor + expect_equal( + names(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)), + c("exp", "obs", "exp_cor") + ) + expect_equal( + dim(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$exp_cor), + c(sdate = 1, member = 3, dataset = 1) + ) + expect_equal( + dim(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$exp), + c(sdate = 3, member = 2, dataset = 1) + ) + expect_equal( + dim(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$obs), + c(sdate = 3, member = 1, dataset = 1) + ) + expect_equal( + c(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$exp_cor), + c(0.3896168, 0.4384543, -0.1302738), + tolerance = 0.0001 + ) + expect_equal( + c(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$exp), + c(-0.1688756, -0.2658420, -0.8049575, -0.3022108, -0.3655258, 0.1237722), + tolerance = 0.0001 + ) + expect_equal( + c(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$obs), + c(-0.1762489, 0.1364694, -1.4581406), + tolerance = 0.0001 + ) }) ############################################## @@ -256,39 +290,34 @@ test_that("3. dat2", { # exp_cor expect_equal( - dim(NAO(exp = NULL, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp_cor), - c(sdate = 1, member = 2) + names(NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)), + c("exp", "obs", "exp_cor") ) expect_equal( - dim(NAO(exp = NULL, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$obs), - c(sdate = 3) + dim(NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp_cor), + c(sdate = 1, member = 2) ) expect_equal( - c((NAO(exp = NULL, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp_cor)), - c(0.2121340, 0.1634516), - tolerance = 0.0001) - expect_equal( - c((NAO(exp = NULL, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$obs)), - c(0.3511294, -0.7196260, -1.5123894), - tolerance = 0.0001 + dim(NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp), + c(sdate = 3, member = 2) ) - expect_equal( - dim(NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)$exp_cor), - c(sdate = 1, member = 2) + dim(NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$obs), + c(sdate = 3) ) expect_equal( - dim(NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)$exp), - c(sdate = 3, member = 2) + c((NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp_cor)), + c(0.2121340, 0.1634516), + tolerance = 0.0001 ) expect_equal( - c(NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)$exp), - c(-0.4225802, -1.3616918, 0.2808729, 0.2723077, -0.7584804, -0.7264514), + c((NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp)), + c(0.01457391, 0.06668166, 0.20193275, -0.20154315, -0.49487925, -0.04181974), tolerance = 0.0001 - ) + ) expect_equal( - c(NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)$exp_cor), - c(-0.22949531, -0.06946422), + c((NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$obs)), + c(0.3511294, -0.7196260, -1.5123894), tolerance = 0.0001 ) -- GitLab From 5e75233a47c0182f498f6b3cc43c8e0b5fe24910 Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 15 Nov 2023 14:57:25 +0100 Subject: [PATCH 11/12] Change T to TRUE --- R/NAO.R | 2 +- man/NAO.Rd | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/NAO.R b/R/NAO.R index 2a0e118..0810f2b 100644 --- a/R/NAO.R +++ b/R/NAO.R @@ -88,7 +88,7 @@ #'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 = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon, obsproj = T) +#'nao <- NAO(exp = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon, obsproj = TRUE) #'# plot the NAO index #' \dontrun{ #'nao$exp <- Reorder(nao$exp, c(2, 1)) diff --git a/man/NAO.Rd b/man/NAO.Rd index 8545e6b..8115fa2 100644 --- a/man/NAO.Rd +++ b/man/NAO.Rd @@ -114,7 +114,7 @@ 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 = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon, obsproj = T) +nao <- NAO(exp = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon, obsproj = TRUE) # plot the NAO index \dontrun{ nao$exp <- Reorder(nao$exp, c(2, 1)) -- GitLab From 8fb2a2e8576e4efc87c926e5070e3384eefc545d Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 15 Nov 2023 15:27:35 +0100 Subject: [PATCH 12/12] Fix example for new NAO dev --- R/PlotBoxWhisker.R | 2 +- R/PlotStereoMap.R | 2 ++ man/PlotBoxWhisker.Rd | 2 +- man/PlotStereoMap.Rd | 2 ++ 4 files changed, 6 insertions(+), 2 deletions(-) diff --git a/R/PlotBoxWhisker.R b/R/PlotBoxWhisker.R index 2ddcec0..9b4c88e 100644 --- a/R/PlotBoxWhisker.R +++ b/R/PlotBoxWhisker.R @@ -90,7 +90,7 @@ #'ano <- Ano_CrossValid(sampleData$mod, sampleData$obs) #'ano_exp <- array(ano$exp, dim = dim(ano$exp)[-2]) #'ano_obs <- array(ano$obs, dim = dim(ano$obs)[-2]) -#'nao <- NAO(ano_exp, ano_obs, sampleData$lat, sampleData$lon) +#'nao <- NAO(exp = ano_exp, obs = ano_obs, lat = sampleData$lat, lon = sampleData$lon) #'# Finally plot the nao index #' \dontrun{ #'nao$exp <- Reorder(nao$exp, c(2, 1)) diff --git a/R/PlotStereoMap.R b/R/PlotStereoMap.R index 4b4fbd2..a2fdb85 100644 --- a/R/PlotStereoMap.R +++ b/R/PlotStereoMap.R @@ -170,8 +170,10 @@ #'data <- matrix(rnorm(100 * 50), 100, 50) #'x <- seq(from = 0, to = 360, length.out = 100) #'y <- seq(from = -90, to = 90, length.out = 50) +#' \dontrun{ #'PlotStereoMap(data, x, y, latlims = c(60, 90), brks = 50, #' toptitle = "This is the title") +#' } #'@import mapproj #'@importFrom grDevices dev.cur dev.new dev.off gray #'@importFrom stats median diff --git a/man/PlotBoxWhisker.Rd b/man/PlotBoxWhisker.Rd index 9c5a3f4..e868160 100644 --- a/man/PlotBoxWhisker.Rd +++ b/man/PlotBoxWhisker.Rd @@ -123,7 +123,7 @@ sampleData$lat[] <- c(20, 80) ano <- Ano_CrossValid(sampleData$mod, sampleData$obs) ano_exp <- array(ano$exp, dim = dim(ano$exp)[-2]) ano_obs <- array(ano$obs, dim = dim(ano$obs)[-2]) -nao <- NAO(ano_exp, ano_obs, sampleData$lat, sampleData$lon) +nao <- NAO(exp = ano_exp, obs = ano_obs, lat = sampleData$lat, lon = sampleData$lon) # Finally plot the nao index \dontrun{ nao$exp <- Reorder(nao$exp, c(2, 1)) diff --git a/man/PlotStereoMap.Rd b/man/PlotStereoMap.Rd index 1b7f166..10e7503 100644 --- a/man/PlotStereoMap.Rd +++ b/man/PlotStereoMap.Rd @@ -283,6 +283,8 @@ compatible with figure layouts if colour bar is disabled. data <- matrix(rnorm(100 * 50), 100, 50) x <- seq(from = 0, to = 360, length.out = 100) y <- seq(from = -90, to = 90, length.out = 50) + \dontrun{ PlotStereoMap(data, x, y, latlims = c(60, 90), brks = 50, toptitle = "This is the title") } +} -- GitLab