Commit 33c5e84e authored by Carlos Delgado Torres's avatar Carlos Delgado Torres
Browse files

removed member_dim parameter. added na.rm parameter. now it uses multiApply

parent 84613b95
Pipeline #5154 failed with stage
in 2 minutes and 39 seconds
......@@ -8,11 +8,10 @@
#'weighted-averaged SST anomalies over 60ºS-60ºN, 0ºE-360ºE (Trenberth &
#'Dennis, 2005; Doblas-Reyes et al., 2013).
#'
#'@param data A numerical array to be used for the index computation with the
#' dimensions: 1) latitude, longitude, start date, forecast month, and member
#' (in case of decadal predictions), 2) latitude, longitude, year, month and
#' member (in case of historical simulations), or 3) latitude, longitude, year
#' and month (in case of observations or reanalyses). This data has to be
#'@param data A numerical array to be used for the index computation with, at least, the
#' dimensions: 1) latitude, longitude, start date and forecast month
#' (in case of decadal predictions), 2) latitude, longitude, year and month
#' (in case of historical simulations or observations). This data has to be
#' provided, at least, over the whole region needed to compute the index.
#'@param data_lats A numeric vector indicating the latitudes of the data.
#'@param data_lons A numeric vector indicating the longitudes of the data.
......@@ -45,7 +44,7 @@
#' anomalies, set it to FALSE. The default value is NULL.\cr
#' In case of parameter 'type' is 'dcpp', 'indices_for_clim' must be relative
#' to the first forecast year, and the climatology is automatically computed
#' over the actual common period for the different forecast years.
#' over the common calendar period for the different forecast years.
#'@param year_dim A character string indicating the name of the year dimension
#' The default value is 'year'. Only used if parameter 'type' is 'hist' or
#' 'obs'.
......@@ -55,13 +54,14 @@
#'@param member_dim A character string indicating the name of the member
#' dimension. The default value is 'member'. Only used if parameter 'type' is
#' 'dcpp' or 'hist'.
#'@param na.rm A logical value indicanting whether to remove NA values. 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 numerical array of the AMV index with the dimensions of:
#' 1) sdate, forecast year, and member (in case of decadal predictions);
#' 2) year and member (in case of historical simulations); or
#' 3) year (in case of observations or reanalyses).
#'@return A numerical array with the AMV index with the same dimensions as data except
#' the lat_dim, lon_dim and fmonth_dim (month_dim) in case of decadal predictions
#' (historical simulations or observations)
#'
#'@examples
#' ## Observations or reanalyses
......@@ -88,7 +88,7 @@
AMV <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'lon',
mask = NULL, monini = 11, fmonth_dim = 'fmonth', sdate_dim = 'sdate',
indices_for_clim = NULL, year_dim = 'year', month_dim = 'month',
member_dim = 'member', ncores = NULL) {
na.rm = TRUE, ncores = NULL) {
## Input Checks
# data
......@@ -209,16 +209,11 @@ AMV <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'lo
stop("Parameter 'month_dim' is not found in 'data' dimension.")
}
}
# member_dim
if (type == 'hist' | type == 'dcpp') {
if (!(is.character(member_dim) & length(member_dim) == 1)) {
stop("Parameter 'member_dim' must be a character string.")
}
if (!member_dim %in% names(dim(data))) {
stop("Parameter 'member_dim' is not found in 'data' dimension.")
}
# na.rm
if (!na.rm %in% c(TRUE,FALSE)) {
stop("Parameter 'na.rm' must be TRUE or FALSE")
}
## Regions for AMV (Doblas-Reyes et al., 2013)
lat_min_1 <- 0; lat_max_1 <- 60
lon_min_1 <- 280; lon_max_1 <- 359.9
......@@ -239,9 +234,17 @@ AMV <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'lo
data <- ClimProjDiags::CombineIndices(indices = list(mean_1,mean_2), weights = NULL,
operation = 'subtract') # (mean_1 - mean_2)
INDEX <- .Indices(data = data, type = type, monini = monini,
indices_for_clim = indices_for_clim, fmonth_dim = fmonth_dim,
sdate_dim = sdate_dim, year_dim = year_dim,
month_dim = month_dim, member_dim = member_dim)
if (type == 'dcpp'){
target_dims <- c(sdate_dim, fmonth_dim)
} else if (type %in% c('hist','obs')){
target_dims <- c(year_dim, month_dim)
}
source('/esarchive/scratch/cdelgado/gitlab/s2dv/R/Utils.R')
INDEX <- multiApply::Apply(data = data, target_dims = target_dims, fun = .Indices,
type = type, monini = monini, indices_for_clim = indices_for_clim,
fmonth_dim = fmonth_dim, sdate_dim = sdate_dim,
year_dim = year_dim, month_dim = month_dim,
na.rm = na.rm, ncores = ncores)$output1
return(INDEX)
}
......@@ -1670,107 +1670,59 @@
}
# to be used in AMV.R, TPI.R, SPOD.R, GSAT.R and GMST.R
.Indices <- function(data, type, monini, indices_for_clim,
fmonth_dim, sdate_dim, year_dim, month_dim, member_dim) {
.Indices <- function(data, type, monini, indices_for_clim,
fmonth_dim, sdate_dim, year_dim, month_dim, na.rm) {
data = drop(data)
if(member_dim %in% names(dim(data))){
if (type == 'dcpp'){
data = s2dv::Reorder(data = data, order = c(sdate_dim,fmonth_dim,member_dim))
} else if (type %in% c('hist','obs')){
data = s2dv::Reorder(data = data, order = c(year_dim,month_dim,member_dim))
}
}
if (type == 'dcpp'){
if (type == 'dcpp') {
data = s2dv::Season(data = data, time_dim = fmonth_dim,
monini = monini, moninf = 1, monsup = 12,
method = mean, na.rm = FALSE)
names(dim(data))[which(names(dim(data))==fmonth_dim)] = 'fyear'
if (member_dim %in% names(dim(data))){
data = s2dv::Reorder(data = data, order = c('fyear',sdate_dim,member_dim))
} else {
data = s2dv::Reorder(data = data, order = c('fyear',sdate_dim))
}
fyear_dim <- 'fyear'
data <- s2dv::Season(data = data, time_dim = fmonth_dim,
monini = monini, moninf = 1, monsup = 12,
method = mean, na.rm = na.rm)
names(dim(data))[which(names(dim(data))==fmonth_dim)] <- fyear_dim
if (is.logical(indices_for_clim)) {
if(!any(indices_for_clim)) {
# indices_for_clim == FALSE -> anomalies are directly given
anom = data
}
if (identical(indices_for_clim, FALSE)) { ## data is already anomalies
} else {
## Different indices_for_clim for each forecast year (same actual years)
anom <- data
n_fyears = as.numeric(dim(data)['fyear'])
n_sdates = as.numeric(dim(data)[sdate_dim])
} else { ## Different indices_for_clim for each forecast year (to use the same calendar years)
if (is.null(indices_for_clim)){
# indices_for_clim == NULL -> anomalies based on the whole (common) period
first_years_for_clim = n_fyears : 1
last_years_for_clim = n_sdates : (n_sdates - n_fyears + 1)
} else {
first_years_for_clim = seq(from = indices_for_clim[1], by = -1, length.out = n_fyears)
last_years_for_clim = seq(from = indices_for_clim[length(indices_for_clim)], by = -1, length.out = n_fyears)
}
anom = array(data = NA, dim = dim(data))
if (member_dim %in% names(dim(data))){
clim = array(data = NA, dim = c(dim(data)['fyear'],dim(data)[member_dim]))
} else {
clim = array(data = NA, dim = c(dim(data)['fyear']))
n_fyears <- as.numeric(dim(data)[fyear_dim])
n_sdates <- as.numeric(dim(data)[sdate_dim])
if (is.null(indices_for_clim)) { ## climatology over the whole (common) period
first_years_for_clim <- n_fyears : 1
last_years_for_clim <- n_sdates : (n_sdates - n_fyears + 1)
} else { ## indices_for_clim specified as a numeric vector
first_years_for_clim <- seq(from = indices_for_clim[1], by = -1, length.out = n_fyears)
last_years_for_clim <- seq(from = indices_for_clim[length(indices_for_clim)], by = -1, length.out = n_fyears)
}
for (i in 1:n_fyears){
if (member_dim %in% names(dim(data))){
for (m in 1:as.numeric(dim(data)[member_dim])){
clim[i,m] = mean(data[i,first_years_for_clim[i]:last_years_for_clim[i],m])
anom[i,,m] = data[i,,m] - clim[i,m]
}
} else {
clim = mean(data[i,first_years_for_clim[i]:last_years_for_clim[i]])
anom[i,] = data[i,] - clim
}
data <- s2dv::Reorder(data = data, order = c(fyear_dim, sdate_dim))
anom <- array(data = NA, dim = dim(data))
clim <- array(data = NA, dim = c(dim(data)[fyear_dim]))
for (i in 1:n_fyears) {
clim <- mean(data[i,first_years_for_clim[i]:last_years_for_clim[i]], na.rm = na.rm)
anom[i,] <- data[i,] - clim
}
}
} else if (type %in% c('obs','hist')){
} else if (type %in% c('obs','hist')) {
data = multiApply::Apply(data = data, target_dims = month_dim, fun = mean)$output1
data <- multiApply::Apply(data = data, target_dims = month_dim, fun = mean, na.rm = na.rm)$output1
if (is.logical(indices_for_clim)) {
if(!any(indices_for_clim)) {
anom = data
}
} else {
if (is.null(indices_for_clim)){
clim = multiApply::Apply(data = data, target_dims = year_dim, fun = mean)$output1
} else {
if (member_dim %in% names(dim(data))){
target_dims = c(year_dim,member_dim)
} else {
target_dims = year_dim
}
clim = multiApply::Apply(data = ClimProjDiags::Subset(x = data, along = year_dim, indices = indices_for_clim),
target_dims = target_dims, fun = mean)$output1
}
anom = multiApply::Apply(data = data, target_dims = year_dim,
fun = function(data,clim){data-clim}, clim = clim)$output1
if (identical(indices_for_clim, FALSE)) { ## data is already anomalies
clim <- 0
} else if (is.null(indices_for_clim)) { ## climatology over the whole period
clim <- multiApply::Apply(data = data, target_dims = year_dim, fun = mean, na.rm = na.rm)$output1
} else { ## indices_for_clim specified as a numeric vector
clim <- multiApply::Apply(data = ClimProjDiags::Subset(x = data, along = year_dim, indices = indices_for_clim),
target_dims = year_dim, fun = mean, na.rm = na.rm)$output1
}
anom <- data - clim
} else {stop('type must be dcpp, hist or obs')}
return(anom)
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment