From b7867b651263b7fdd2557886c1c20efecf4658e3 Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Fri, 19 Jun 2020 16:26:47 +0200 Subject: [PATCH 01/11] added the AMV.R, GMST.R, GSAT.R, SPOD.R, TPI.R files and the .Indices function in Utils.R --- R/AMV.R | 140 +++++++++++++++++++++++++++++++++++++++++++++ R/GMST.R | 166 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ R/GSAT.R | 125 ++++++++++++++++++++++++++++++++++++++++ R/SPOD.R | 139 +++++++++++++++++++++++++++++++++++++++++++++ R/TPI.R | 145 +++++++++++++++++++++++++++++++++++++++++++++++ R/Utils.R | 99 ++++++++++++++++++++++++++++++++ 6 files changed, 814 insertions(+) create mode 100644 R/AMV.R create mode 100644 R/GMST.R create mode 100644 R/GSAT.R create mode 100644 R/SPOD.R create mode 100644 R/TPI.R diff --git a/R/AMV.R b/R/AMV.R new file mode 100644 index 0000000..f5b70e5 --- /dev/null +++ b/R/AMV.R @@ -0,0 +1,140 @@ +#'Compute the Atlantic Multidecadal Variability (AMV) index +#' +#'@description The Atlantic Multidecadal Variability (AMV), also known as Atlantic Multidecadal Oscillation (AMO), +#'is a model of natural variability of the sea surface temperatures (SST) over the North Atlantic Ocean +#'on multi-decadal time scales. The AMV index is computed as the difference of weighted-averaged SST anomalies over +#'the North Atlantic region (0ºN-60ºN, 280ºE-360ºE) and the weighted-averaged SST anomalies over 60ºS-60ºN, 0ºE-360ºE +#'(Trenberth & Dennis, 2005; Doblas-Reyes et al., 2013). +#'@author Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} +#'@author Roberto Bilbao, \email{roberto.bilbao@bsc.es} +#'@author Núria Pérez-Zanón, \email{nuria.perez@bsc.es} +#' +#'@param data Data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), +#' with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). +#' This data has to be provided, at least, over the whole region needed to compute the index. +#'@param data_lats An array with the latitudes of the data. +#'@param data_lons An array with the longitudes of the data. +#'@param type A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses). +#'@param lat_dim A string with the name of the latitude dimension ('lat' by default). +#'@param lon_dim A string with the name of the longitude dimension ('lon' by default). +#'@param mask An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a +#' combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the +#' observational dataset for a fair comparison between the forecast system and the reference dataset. +#'@param monini Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'. +#'@param fmonth_dim A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'. +#'@param sdate_dim A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'. +#'@param indices_for_clim Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). +#' In case of type='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. +#'@param year_dim A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'. +#'@param month_dim A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'. +#'@param member_dim A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'. +#' +#'@return The AMV index as funcion of the sdate, forecast year, and member (in case of decadal predictions); +#' as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +#' +#'@examples +#' ## Observations or reanalyses +#' obs = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +#' lat = seq(-90, 90, 10) +#' lon = seq(0, 360, 10) +#' index_obs = AMV(data = obs, data_lats = lat, data_lons = lon, type = 'obs') +#' +#' ## Historical simulations +#' hist = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +#' lat = seq(-90, 90, 10) +#' lon = seq(0, 360, 10) +#' index_hist = AMV(data = hist, data_lats = lat, data_lons = lon, type = 'hist') +#' +#' ## Decadal predictions +#' dcpp = array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +#' lat = seq(-90, 90, 10) +#' lon = seq(0, 360, 10) +#' index_dcpp = AMV(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) +#' +#'@import ClimProjDiags +#'@import multiApply +#'@import s2dv +#'@import startR +#'@export +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'){ + + library(ClimProjDiags) + library(multiApply) + library(s2dv) + library(startR) + source('/esarchive/scratch/cdelgado/gitlab/c3s-scripts/R/Indices.R') + + ## Checkings + if (!is.array(data)){ + stop('data must be an array') + } + if(!class(data_lats)=='numeric'){ + stop('data_lats must be a numeric vector') + } + if (!class(data_lons)=='numeric'){ + stop('data_lons must be a numeric vector or NULL') + } + if (!type %in% c('dcpp','hist','obs')){ + stop("type must be 'dcpp', 'hist' or 'obs'") + } + if (!is.character(lat_dim)){ + stop('lat_dim must be a string') + } + if (!is.character(lon_dim)){ + stop('lon_dim must be a string') + } + if (!monini){ + stop("monini must be an integer from 1 to 12") + } + if (!is.character(fmonth_dim)){ + stop('fmonth_dim must be a string') + } + if (!is.character(sdate_dim)){ + stop('sdate_dim must be a string') + } + if (!is.null(indices_for_clim) & !class(indices_for_clim) %in% c('numeric','integer') & !is.null(indices_for_clim) & !isFALSE(indices_for_clim)){ + stop("indices_for_clim must be a numeric vector, NULL to compute the anomalies based on the whole period, or FALSE if data are already anomalies") + } + if (!is.character(year_dim)){ + stop('year_dim must be a string') + } + if (!is.character(month_dim)){ + stop('month_dim must be a string') + } + if (!is.null(mask)){ + if (is.array(mask) & identical(names(dim(mask)),c(lat_dim,lon_dim)) & identical(as.integer(dim(mask)),c(length(data_lats),length(data_lons)))){ + ## To mask those grid point that are missing in the observations + mask <- s2dv::Reorder(data = mask, order = c(lat_dim,lon_dim)) + fun_mask <- function(data, mask){ + data[mask == 0] <- NA + return(data) + } + data <- multiApply::Apply(data = data, target_dims = c(lat_dim,lon_dim), fun = fun_mask, mask = mask)$output1 + } else { + stop('mask must be NULL (no mask) or an array with c(lat_dim,lon_dim) dimensions and 0 in those grid points that have to be masked') + } + } + + ## 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 + lat_min_2 <- 60; lat_max_2 <- -60 + lon_min_2 <- 0; lon_max_2 <- 359.9 + regions <- NULL + regions$reg1 <- c(lon_min_1, lon_max_1, lat_min_1, lat_max_1) + regions$reg2 <- c(lon_min_2, lon_max_2, lat_min_2, lat_max_2) + + mean_1 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = regions$reg1, + londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) + mean_2 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = regions$reg2, + londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) + 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) + return(INDEX) +} \ No newline at end of file diff --git a/R/GMST.R b/R/GMST.R new file mode 100644 index 0000000..81ebb7d --- /dev/null +++ b/R/GMST.R @@ -0,0 +1,166 @@ +#'Compute the Global Mean Surface Temperature (GMST) anomalies +#' +#'@description The Global Mean Surface Temperature (GMST) anomalies are computed as the weighted-averaged +#'surface air temperature anomalies over land and sea surface temperature anomalies over the ocean. +#'@author Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} +#'@author Roberto Bilbao, \email{roberto.bilbao@bsc.es} +#'@author Núria Pérez-Zanón, \email{nuria.perez@bsc.es} +#' +#'@param data_tas Surface air temperature data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), +#' with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). +#' This data has to be provided, at least, over the whole region needed to compute the index. The dimensions must be identical to those of data_tos. +#'@param data_tos Sea surface temperature data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), +#' with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). +#' This data has to be provided, at least, over the whole region needed to compute the index. The dimensions must be identical to those of data_tas. +#'@param data_lats An array with the latitudes of the data. +#'@param data_lons An array with the longitudes of the data. +#'@param mask_sea_land # mask for blending the surface air temperature over land and the sea surface temperature over the ocean. It must have lat_dim and lon_dim dimensions, and their lengths +#' must be equal to the length of data_lats and data_lons. +#'@param sea_value Value of the sea grid points in the mask. +#'@param type A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses). +#'@param lat_dim A string with the name of the latitude dimension ('lat' by default). +#'@param lon_dim A string with the name of the longitude dimension ('lon' by default). +#'@param mask An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a +#' combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the +#' observational dataset for a fair comparison between the forecast system and the reference dataset. +#'@param monini Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'. +#'@param fmonth_dim A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'. +#'@param sdate_dim A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'. +#'@param indices_for_clim Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). +#' In case of type='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. +#'@param year_dim A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'. +#'@param month_dim A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'. +#'@param member_dim A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'. +#' +#'@return The GMST anomalies as function of the sdate, forecast year, and member (in case of decadal predictions); +#' as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +#' +#'@examples +#' ## Observations or reanalyses +#' obs_tas = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +#' obs_tos = array(2:101, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +#' mask_sea_land = array(c(1,0,1), dim = c(lat = 19, lon = 37)) +#' sea_value = 1 +#' lat = seq(-90, 90, 10) +#' lon = seq(0, 360, 10) +#' index_obs = GMST(data_tas = obs_tas, data_tos = obs_tos, data_lats = lat, data_lons = lon, type = 'obs', mask_sea_land = mask_sea_land, sea_value = sea_value) +#' +#' ## Historical simulations +#' hist_tas = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +#' hist_tos = array(2:101, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +#' mask_sea_land = array(c(1,0,1), dim = c(lat = 19, lon = 37)) +#' sea_value = 1 +#' lat = seq(-90, 90, 10) +#' lon = seq(0, 360, 10) +#' index_hist = GMST(data_tas = hist_tas, data_tos = hist_tos, data_lats = lat, data_lons = lon, type = 'hist', mask_sea_land = mask_sea_land, sea_value = sea_value) +#' +#' ## Decadal predictions +#' dcpp_tas = array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +#' dcpp_tos = array(2:101, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +#' mask_sea_land = array(c(1,0,1), dim = c(lat = 19, lon = 37)) +#' sea_value = 1 +#' lat = seq(-90, 90, 10) +#' lon = seq(0, 360, 10) +#' index_dcpp = GMST(data_tas = dcpp_tas, data_tos = dcpp_tos, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1, mask_sea_land = mask_sea_land, sea_value = sea_value) +#' +#'@import ClimProjDiags +#'@import multiApply +#'@import s2dv +#'@import startR +#'@export +GMST <- function(data_tas, data_tos, data_lats, data_lons, mask_sea_land, sea_value, + type, mask = NULL, lat_dim = 'lat', lon_dim = 'lon', monini = 11, + fmonth_dim = 'fmonth', sdate_dim = 'sdate', indices_for_clim = NULL, + year_dim = 'year', month_dim = 'month', member_dim = 'member'){ + + library(ClimProjDiags) + library(multiApply) + library(s2dv) + library(startR) + source('/esarchive/scratch/cdelgado/gitlab/c3s-scripts/R/Indices.R') + + ## Checkings + if (!is.array(data_tas)){ + stop('data_tas must be an array') + } + if (!is.array(data_tos)){ + stop('data_tos must be an array') + } + if (!identical(dim(data_tas),dim(data_tos))){ + stop('The dimension of data_tas and data_tos must be identical') + } + if(!class(data_lats)=='numeric'){ + stop('data_lats must be a numeric vector') + } + if (!class(data_lons)=='numeric'){ + stop('data_lons must be a numeric vector or NULL') + } + if (!type %in% c('dcpp','hist','obs')){ + stop("type must be 'dcpp', 'hist' or 'obs'") + } + if (!is.character(lat_dim)){ + stop('lat_dim must be a string') + } + if (!is.character(lon_dim)){ + stop('lon_dim must be a string') + } + if (!monini){ + stop("monini must be an integer from 1 to 12") + } + if (!is.character(fmonth_dim)){ + stop('fmonth_dim must be a string') + } + if (!is.character(sdate_dim)){ + stop('sdate_dim must be a string') + } + if (!is.null(indices_for_clim) & !class(indices_for_clim) %in% c('numeric','integer') & !is.null(indices_for_clim) & !isFALSE(indices_for_clim)){ + stop("indices_for_clim must be a numeric vector, NULL to compute the anomalies based on the whole period, or FALSE if data are already anomalies") + } + if (!is.character(year_dim)){ + stop('year_dim must be a string') + } + if (!is.character(month_dim)){ + stop('month_dim must be a string') + } + if (!is.array(mask_sea_land)){ + stop('mask_sea_land must be an array with c(lat_dim,lon_dim) dimensions') + } else if (!identical(names(dim(mask_sea_land)),c(lat_dim,lon_dim))){ + stop('mask_sea_land must be an array with c(lat_dim,lon_dim) dimensions') + } else if (!identical(as.integer(dim(mask_sea_land)),c(length(data_lats),length(data_lons)))){ + stop('mask_sea_land dimension must be equal to the data_lats and data_lons lengths') + } + ## combination of tas and tos (data) + mask_tas_tos <- function(data_tas, data_tos, mask_sea_land, sea_value) { + data <- data_tas + data[mask_sea_land==sea_value] <- data_tos[mask_sea_land==sea_value] + return(data) + } + mask_sea_land <- s2dv::Reorder(data = mask_sea_land, order = c(lat_dim,lon_dim)) + data <- multiApply::Apply(data = list(data_tas, data_tos), target_dims = c(lat_dim,lon_dim), + fun = mask_tas_tos, mask_sea_land = mask_sea_land, sea_value = sea_value)$output1 + data <- drop(data) + rm(data_tas, data_tos) + + if (!is.null(mask)){ + if (is.array(mask) & identical(names(dim(mask)),c(lat_dim,lon_dim)) & identical(as.integer(dim(mask)),c(length(data_lats),length(data_lons)))){ + ## To mask those grid point that are missing in the observations + mask <- s2dv::Reorder(data = mask, order = c(lat_dim,lon_dim)) + fun_mask <- function(data, mask){ + data[mask == 0] <- NA + return(data) + } + data <- multiApply::Apply(data = data, target_dims = c(lat_dim,lon_dim), fun = fun_mask, mask = mask)$output1 + } else { + stop('mask must be NULL (no mask) or an array with c(lat_dim,lon_dim) dimensions and 0 in those grid points that have to be masked') + } + } + + data <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = NULL, + londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) + + 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) + return(INDEX) +} \ No newline at end of file diff --git a/R/GSAT.R b/R/GSAT.R new file mode 100644 index 0000000..e5bce2a --- /dev/null +++ b/R/GSAT.R @@ -0,0 +1,125 @@ +#'Compute the Global Surface Air Temperature (GSAT) anomalies +#' +#'@description The Global Surface Air Temperature (GSAT) anomalies are computed as the +#'weighted-averaged surface air temperature anomalies over the global region. +#'@author Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} +#'@author Núria Pérez-Zanón, \email{nuria.perez@bsc.es} +#'@author Roberto Bilbao, \email{roberto.bilbao@bsc.es} +#' +#'@param data Data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), +#' with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). +#' This data has to be provided, at least, over the whole region needed to compute the index. +#'@param data_lats An array with the latitudes of the data. +#'@param data_lons An array with the longitudes of the data. +#'@param type A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses). +#'@param lat_dim A string with the name of the latitude dimension ('lat' by default). +#'@param lon_dim A string with the name of the longitude dimension ('lon' by default). +#'@param mask An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a +#' combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the +#' observational dataset for a fair comparison between the forecast system and the reference dataset. +#'@param monini Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'. +#'@param fmonth_dim A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'. +#'@param sdate_dim A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'. +#'@param indices_for_clim Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). +#' In case of type='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. +#'@param year_dim A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'. +#'@param month_dim A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'. +#'@param member_dim A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'. +#' +#'@return The GSAT anomalies as function of the sdate, forecast year, and member (in case of decadal predictions); +#' as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +#' +#'@examples +#' ## Observations or reanalyses +#' obs = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +#' lat = seq(-90, 90, 10) +#' lon = seq(0, 360, 10) +#' index_obs = GSAT(data = obs, data_lats = lat, data_lons = lon, type = 'obs') +#' +#' ## Historical simulations +#' hist = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +#' lat = seq(-90, 90, 10) +#' lon = seq(0, 360, 10) +#' index_hist = GSAT(data = hist, data_lats = lat, data_lons = lon, type = 'hist') +#' +#' ## Decadal predictions +#' dcpp = array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +#' lat = seq(-90, 90, 10) +#' lon = seq(0, 360, 10) +#' index_dcpp = GSAT(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) +#' +#'@import ClimProjDiags +#'@import multiApply +#'@import s2dv +#'@import startR +#'@export +GSAT <- 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'){ + + library(ClimProjDiags) + library(multiApply) + library(s2dv) + library(startR) + source('/esarchive/scratch/cdelgado/gitlab/c3s-scripts/R/Indices.R') + + ## Checkings + if (!is.array(data)){ + stop('data must be an array') + } + if(!class(data_lats)=='numeric'){ + stop('data_lats must be a numeric vector') + } + if (!class(data_lons)=='numeric'){ + stop('data_lons must be a numeric vector or NULL') + } + if (!type %in% c('dcpp','hist','obs')){ + stop("type must be 'dcpp', 'hist' or 'obs'") + } + if (!is.character(lat_dim)){ + stop('lat_dim must be a string') + } + if (!is.character(lon_dim)){ + stop('lon_dim must be a string') + } + if (!monini){ + stop("monini must be an integer from 1 to 12") + } + if (!is.character(fmonth_dim)){ + stop('fmonth_dim must be a string') + } + if (!is.character(sdate_dim)){ + stop('sdate_dim must be a string') + } + if (!is.null(indices_for_clim) & !class(indices_for_clim) %in% c('numeric','integer') & !is.null(indices_for_clim) & !isFALSE(indices_for_clim)){ + stop("indices_for_clim must be a numeric vector, NULL to compute the anomalies based on the whole period, or FALSE if data are already anomalies") + } + if (!is.character(year_dim)){ + stop('year_dim must be a string') + } + if (!is.character(month_dim)){ + stop('month_dim must be a string') + } + if (!is.null(mask)){ + if (is.array(mask) & identical(names(dim(mask)),c(lat_dim,lon_dim)) & identical(as.integer(dim(mask)),c(length(data_lats),length(data_lons)))){ + ## To mask those grid point that are missing in the observations + mask <- s2dv::Reorder(data = mask, order = c(lat_dim,lon_dim)) + fun_mask <- function(data, mask){ + data[mask == 0] <- NA + return(data) + } + data <- multiApply::Apply(data = data, target_dims = c(lat_dim,lon_dim), fun = fun_mask, mask = mask)$output1 + } else { + stop('mask must be NULL (no mask) or an array with c(lat_dim,lon_dim) dimensions and 0 in those grid points that have to be masked') + } + } + + data <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = NULL, + londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) + + 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) + return(INDEX) +} \ No newline at end of file diff --git a/R/SPOD.R b/R/SPOD.R new file mode 100644 index 0000000..0e4cc8d --- /dev/null +++ b/R/SPOD.R @@ -0,0 +1,139 @@ +#'Compute the South Pacific Ocean Dipole (SPOD) index +#' +#'@description The South Pacific Ocean Dipole (SPOD) index is related to the El Nino-Southern Oscillation (ENSO) and the Inderdecadal Pacific Oscillation (IPO). +#'The SPOD index is computed as the difference of weighted-averaged SST anomalies over 20ºS-48ºS, 165ºE-190ºE (NW pole) and the weighted-averaged SST anomalies over +#'44ºS-65ºS, 220ºE-260ºE (SE pole) (Saurral et al., 2020). +#'@author Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} +#'@author Núria Pérez-Zanón, \email{nuria.perez@bsc.es} +#'@author Roberto Bilbao, \email{roberto.bilbao@bsc.es} +#' +#'@param data Data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), +#' with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). +#' This data has to be provided, at least, over the whole region needed to compute the index. +#'@param data_lats An array with the latitudes of the data. +#'@param data_lons An array with the longitudes of the data. +#'@param type A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses). +#'@param lat_dim A string with the name of the latitude dimension ('lat' by default). +#'@param lon_dim A string with the name of the longitude dimension ('lon' by default). +#'@param mask An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a +#' combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the +#' observational dataset for a fair comparison between the forecast system and the reference dataset. +#'@param monini Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'. +#'@param fmonth_dim A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'. +#'@param sdate_dim A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'. +#'@param indices_for_clim Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). +#' In case of type='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. +#'@param year_dim A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'. +#'@param month_dim A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'. +#'@param member_dim A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'. +#' +#'@return The SPOD index as function of the sdate, forecast year, and member (in case of decadal predictions); +#' as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +#' +#'@examples +#' ## Observations or reanalyses +#' obs = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +#' lat = seq(-90, 90, 10) +#' lon = seq(0, 360, 10) +#' index_obs = SPOD(data = obs, data_lats = lat, data_lons = lon, type = 'obs') +#' +#' ## Historical simulations +#' hist = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +#' lat = seq(-90, 90, 10) +#' lon = seq(0, 360, 10) +#' index_hist = SPOD(data = hist, data_lats = lat, data_lons = lon, type = 'hist') +#' +#' ## Decadal predictions +#' dcpp = array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +#' lat = seq(-90, 90, 10) +#' lon = seq(0, 360, 10) +#' index_dcpp = SPOD(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) +#' +#'@import ClimProjDiags +#'@import multiApply +#'@import s2dv +#'@import startR +#'@export +SPOD <- 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'){ + + library(ClimProjDiags) + library(multiApply) + library(s2dv) + library(startR) + source('/esarchive/scratch/cdelgado/gitlab/c3s-scripts/R/Indices.R') + + ## Checkings + if (!is.array(data)){ + stop('data must be an array') + } + if(!class(data_lats)=='numeric'){ + stop('data_lats must be a numeric vector') + } + if (!class(data_lons)=='numeric'){ + stop('data_lons must be a numeric vector or NULL') + } + if (!type %in% c('dcpp','hist','obs')){ + stop("type must be 'dcpp', 'hist' or 'obs'") + } + if (!is.character(lat_dim)){ + stop('lat_dim must be a string') + } + if (!is.character(lon_dim)){ + stop('lon_dim must be a string') + } + if (!monini){ + stop("monini must be an integer from 1 to 12") + } + if (!is.character(fmonth_dim)){ + stop('fmonth_dim must be a string') + } + if (!is.character(sdate_dim)){ + stop('sdate_dim must be a string') + } + if (!is.null(indices_for_clim) & !class(indices_for_clim) %in% c('numeric','integer') & !is.null(indices_for_clim) & !isFALSE(indices_for_clim)){ + stop("indices_for_clim must be a numeric vector, NULL to compute the anomalies based on the whole period, or FALSE if data are already anomalies") + } + if (!is.character(year_dim)){ + stop('year_dim must be a string') + } + if (!is.character(month_dim)){ + stop('month_dim must be a string') + } + if (!is.null(mask)){ + if (is.array(mask) & identical(names(dim(mask)),c(lat_dim,lon_dim)) & identical(as.integer(dim(mask)),c(length(data_lats),length(data_lons)))){ + ## To mask those grid point that are missing in the observations + mask <- s2dv::Reorder(data = mask, order = c(lat_dim,lon_dim)) + fun_mask <- function(data, mask){ + data[mask == 0] <- NA + return(data) + } + data <- multiApply::Apply(data = data, target_dims = c(lat_dim,lon_dim), fun = fun_mask, mask = mask)$output1 + } else { + stop('mask must be NULL (no mask) or an array with c(lat_dim,lon_dim) dimensions and 0 in those grid points that have to be masked') + } + } + + ## Regions for IPO_SPOD (Saurral et al., 2020) + lat_min_1 <- -48; lat_max_1 = -20 + lon_min_1 <- 165; lon_max_1 = 190 + lat_min_2 <- -65; lat_max_2 = -44 + lon_min_2 <- 220; lon_max_2 = 260 + regions = NULL + regions$reg1 <- c(lon_min_1, lon_max_1, lat_min_1, lat_max_1) + regions$reg2 <- c(lon_min_2, lon_max_2, lat_min_2, lat_max_2) + + mean_1 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = regions$reg1, + londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) + mean_2 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = regions$reg2, + londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) + + 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) + return(INDEX) +} \ No newline at end of file diff --git a/R/TPI.R b/R/TPI.R new file mode 100644 index 0000000..d0e1efe --- /dev/null +++ b/R/TPI.R @@ -0,0 +1,145 @@ +#'Compute the Tripole Index (TPI) for the Interdecadal Pacific Oscillation (IPO) +#' +#'@description The Tripole Index (TPI) for the Interdecadal Pacific Oscillation (IPO) is computed as the difference of weighted-averaged SST anomalies over +#'10ºS-10ºN, 170ºE-270ºE minus the mean of the weighted-averaged SST anomalies over 25ºN-45ºN, 140ºE-215ºE and 50ºS-15ºS, 150ºE-200ºE (Henley et al., 2015). +#'@author Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} +#'@author Núria Pérez-Zanón, \email{nuria.perez@bsc.es} +#'@author Roberto Bilbao, \email{roberto.bilbao@bsc.es} +#' +#'@param data Data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), +#' with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). +#' This data has to be provided, at least, over the whole region needed to compute the index. +#'@param data_lats An array with the latitudes of the data. +#'@param data_lons An array with the longitudes of the data. +#'@param type A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses). +#'@param lat_dim A string with the name of the latitude dimension ('lat' by default). +#'@param lon_dim A string with the name of the longitude dimension ('lon' by default). +#'@param mask An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a +#' combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the +#' observational dataset for a fair comparison between the forecast system and the reference dataset. +#'@param monini Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'. +#'@param fmonth_dim A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'. +#'@param sdate_dim A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'. +#'@param indices_for_clim Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). +#' In case of type='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. +#'@param year_dim A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'. +#'@param month_dim A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'. +#'@param member_dim A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'. +#' +#'@return The TPI as function of the sdate, forecast year, and member (in case of decadal predictions); +#' as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +#' +#'@examples +#' ## Observations or reanalyses +#' obs = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +#' lat = seq(-90, 90, 10) +#' lon = seq(0, 360, 10) +#' index_obs = TPI(data = obs, data_lats = lat, data_lons = lon, type = 'obs') +#' +#' ## Historical simulations +#' hist = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +#' lat = seq(-90, 90, 10) +#' lon = seq(0, 360, 10) +#' index_hist = TPI(data = hist, data_lats = lat, data_lons = lon, type = 'hist') +#' +#' ## Decadal predictions +#' dcpp = array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +#' lat = seq(-90, 90, 10) +#' lon = seq(0, 360, 10) +#' index_dcpp = TPI(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) +#' +#'@import ClimProjDiags +#'@import multiApply +#'@import s2dv +#'@import startR +#'@export +TPI <- function(data, data_lats, data_lons, index, 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'){ + + library(ClimProjDiags) + library(multiApply) + library(s2dv) + library(startR) + source('/esarchive/scratch/cdelgado/gitlab/c3s-scripts/R/Indices.R') + + ## Checkings + if (!is.array(data)){ + stop('data must be an array') + } + if(!class(data_lats)=='numeric'){ + stop('data_lats must be a numeric vector') + } + if (!class(data_lons)=='numeric'){ + stop('data_lons must be a numeric vector or NULL') + } + if (!type %in% c('dcpp','hist','obs')){ + stop("type must be 'dcpp', 'hist' or 'obs'") + } + if (!is.character(lat_dim)){ + stop('lat_dim must be a string') + } + if (!is.character(lon_dim)){ + stop('lon_dim must be a string') + } + if (!monini){ + stop("monini must be an integer from 1 to 12") + } + if (!is.character(fmonth_dim)){ + stop('fmonth_dim must be a string') + } + if (!is.character(sdate_dim)){ + stop('sdate_dim must be a string') + } + if (!is.null(indices_for_clim) & !class(indices_for_clim) %in% c('numeric','integer') & !is.null(indices_for_clim) & !isFALSE(indices_for_clim)){ + stop("indices_for_clim must be a numeric vector, NULL to compute the anomalies based on the whole period, or FALSE if data are already anomalies") + } + if (!is.character(year_dim)){ + stop('year_dim must be a string') + } + if (!is.character(month_dim)){ + stop('month_dim must be a string') + } + if (!is.null(mask)){ + if (is.array(mask) & identical(names(dim(mask)),c(lat_dim,lon_dim)) & identical(as.integer(dim(mask)),c(length(data_lats),length(data_lons)))){ + ## To mask those grid point that are missing in the observations + mask <- s2dv::Reorder(data = mask, order = c(lat_dim,lon_dim)) + fun_mask <- function(data, mask){ + data[mask == 0] <- NA + return(data) + } + data <- multiApply::Apply(data = data, target_dims = c(lat_dim,lon_dim), fun = fun_mask, mask = mask)$output1 + } else { + stop('mask must be NULL (no mask) or an array with c(lat_dim,lon_dim) dimensions and 0 in those grid points that have to be masked') + } + } + + # Regions for IPO_TPI (psl.noaa.gov/data/timeseries/IPOTPI) + lat_min_1 <- 25; lat_max_1 <- 45 + lon_min_1 <- 140; lon_max_1 <- 215 + lat_min_2 <- -10; lat_max_2 <- 10 + lon_min_2 <- 170; lon_max_2 <- 270 + lat_min_3 <- -50; lat_max_3 <- -15 + lon_min_3 <- 150; lon_max_3 <- 200 + regions = NULL + regions$reg1 <- c(lon_min_1, lon_max_1, lat_min_1, lat_max_1) + regions$reg2 <- c(lon_min_2, lon_max_2, lat_min_2, lat_max_2) + regions$reg3 <- c(lon_min_3, lon_max_3, lat_min_3, lat_max_3) + + mean_1 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = regions$reg1, + londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) + mean_2 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = regions$reg2, + londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) + mean_3 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = regions$reg3, + londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) + + mean_1_3 <- ClimProjDiags::CombineIndices(indices = list(mean_1,mean_3), weights = NULL, operation = 'mean') + + data <- ClimProjDiags::CombineIndices(indices = list(mean_2, mean_1_3), weights = NULL, operation = 'subtract') # mean_2 - ((mean_1 + mean_3)/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) + return(INDEX) +} \ No newline at end of file diff --git a/R/Utils.R b/R/Utils.R index 86b2644..64351c7 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -1660,3 +1660,102 @@ return(output) } +# 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) { + + 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'){ + + 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)) + } + + if (isFALSE(indices_for_clim)){ + + # indices_for_clim == FALSE -> anomalies are directly given + anom = data + + } else { + + ## Different indices_for_clim for each forecast year (same actual years) + + n_fyears = as.numeric(dim(data)['fyear']) + n_sdates = as.numeric(dim(data)[sdate_dim]) + + 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)) + clim = array(data = NA, dim = c(dim(data)['fyear'],dim(data)[member_dim])) + 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 + } + } + } + + } else if (type %in% c('obs','hist')){ + + data = multiApply::Apply(data = data, target_dims = month_dim, fun = mean)$output1 + + if (isFALSE(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 = startR::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 + } + + } else {stop('type must be dcpp, hist or obs')} + + return(anom) +} + -- GitLab From fa74ab99a01a0df094fb99a181469c80a32c6bcf Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Fri, 19 Jun 2020 16:47:54 +0200 Subject: [PATCH 02/11] converted Indices into .Indices --- R/AMV.R | 14 ++++---------- R/GMST.R | 14 ++++---------- R/GSAT.R | 14 ++++---------- R/SPOD.R | 14 ++++---------- R/TPI.R | 14 ++++---------- R/Utils.R | 2 +- 6 files changed, 21 insertions(+), 51 deletions(-) diff --git a/R/AMV.R b/R/AMV.R index f5b70e5..c8f718b 100644 --- a/R/AMV.R +++ b/R/AMV.R @@ -60,12 +60,6 @@ AMV <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'lo mask = NULL, monini = 11, fmonth_dim = 'fmonth', sdate_dim = 'sdate', indices_for_clim = NULL, year_dim = 'year', month_dim = 'month', member_dim = 'member'){ - library(ClimProjDiags) - library(multiApply) - library(s2dv) - library(startR) - source('/esarchive/scratch/cdelgado/gitlab/c3s-scripts/R/Indices.R') - ## Checkings if (!is.array(data)){ stop('data must be an array') @@ -132,9 +126,9 @@ AMV <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'lo londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) 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) + 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) return(INDEX) } \ No newline at end of file diff --git a/R/GMST.R b/R/GMST.R index 81ebb7d..6ea8216 100644 --- a/R/GMST.R +++ b/R/GMST.R @@ -73,12 +73,6 @@ GMST <- function(data_tas, data_tos, data_lats, data_lons, mask_sea_land, sea_va fmonth_dim = 'fmonth', sdate_dim = 'sdate', indices_for_clim = NULL, year_dim = 'year', month_dim = 'month', member_dim = 'member'){ - library(ClimProjDiags) - library(multiApply) - library(s2dv) - library(startR) - source('/esarchive/scratch/cdelgado/gitlab/c3s-scripts/R/Indices.R') - ## Checkings if (!is.array(data_tas)){ stop('data_tas must be an array') @@ -158,9 +152,9 @@ GMST <- function(data_tas, data_tos, data_lats, data_lons, mask_sea_land, sea_va data <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = NULL, londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) - 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) + 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) return(INDEX) } \ No newline at end of file diff --git a/R/GSAT.R b/R/GSAT.R index e5bce2a..8f8d4ab 100644 --- a/R/GSAT.R +++ b/R/GSAT.R @@ -57,12 +57,6 @@ GSAT <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'l mask = NULL, monini = 11, fmonth_dim = 'fmonth', sdate_dim = 'sdate', indices_for_clim = NULL, year_dim = 'year', month_dim = 'month', member_dim = 'member'){ - library(ClimProjDiags) - library(multiApply) - library(s2dv) - library(startR) - source('/esarchive/scratch/cdelgado/gitlab/c3s-scripts/R/Indices.R') - ## Checkings if (!is.array(data)){ stop('data must be an array') @@ -117,9 +111,9 @@ GSAT <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'l data <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = NULL, londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) - 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) + 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) return(INDEX) } \ No newline at end of file diff --git a/R/SPOD.R b/R/SPOD.R index 0e4cc8d..afc5a73 100644 --- a/R/SPOD.R +++ b/R/SPOD.R @@ -58,12 +58,6 @@ SPOD <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'l mask = NULL, monini = 11, fmonth_dim = 'fmonth', sdate_dim = 'sdate', indices_for_clim = NULL, year_dim = 'year', month_dim = 'month', member_dim = 'member'){ - library(ClimProjDiags) - library(multiApply) - library(s2dv) - library(startR) - source('/esarchive/scratch/cdelgado/gitlab/c3s-scripts/R/Indices.R') - ## Checkings if (!is.array(data)){ stop('data must be an array') @@ -131,9 +125,9 @@ SPOD <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'l 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) + 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) return(INDEX) } \ No newline at end of file diff --git a/R/TPI.R b/R/TPI.R index d0e1efe..4e6537e 100644 --- a/R/TPI.R +++ b/R/TPI.R @@ -57,12 +57,6 @@ TPI <- function(data, data_lats, data_lons, index, type, lat_dim = 'lat', lon_di mask = NULL, monini = 11, fmonth_dim = 'fmonth', sdate_dim = 'sdate', indices_for_clim = NULL, year_dim = 'year', month_dim = 'month', member_dim = 'member'){ - library(ClimProjDiags) - library(multiApply) - library(s2dv) - library(startR) - source('/esarchive/scratch/cdelgado/gitlab/c3s-scripts/R/Indices.R') - ## Checkings if (!is.array(data)){ stop('data must be an array') @@ -137,9 +131,9 @@ TPI <- function(data, data_lats, data_lons, index, type, lat_dim = 'lat', lon_di data <- ClimProjDiags::CombineIndices(indices = list(mean_2, mean_1_3), weights = NULL, operation = 'subtract') # mean_2 - ((mean_1 + mean_3)/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) + 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) return(INDEX) } \ No newline at end of file diff --git a/R/Utils.R b/R/Utils.R index 64351c7..2b5569b 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -1661,7 +1661,7 @@ } # to be used in AMV.R, TPI.R, SPOD.R, GSAT.R and GMST.R -Indices <- function(data, type, monini, indices_for_clim, +.Indices <- function(data, type, monini, indices_for_clim, fmonth_dim, sdate_dim, year_dim, month_dim, member_dim) { data = drop(data) -- GitLab From d74e8e9978a8455ac3465ac66a887ecbca13482d Mon Sep 17 00:00:00 2001 From: cdelgado Date: Tue, 23 Jun 2020 18:15:17 +0200 Subject: [PATCH 03/11] Update AMV.R --- R/AMV.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/AMV.R b/R/AMV.R index c8f718b..223edd2 100644 --- a/R/AMV.R +++ b/R/AMV.R @@ -114,7 +114,7 @@ AMV <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'lo ## 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 - lat_min_2 <- 60; lat_max_2 <- -60 + lat_min_2 <- -60; lat_max_2 <- 60 lon_min_2 <- 0; lon_max_2 <- 359.9 regions <- NULL regions$reg1 <- c(lon_min_1, lon_max_1, lat_min_1, lat_max_1) -- GitLab From 8a509f39108fe89c5668301a3431515bd63153d1 Mon Sep 17 00:00:00 2001 From: cdelgado Date: Wed, 1 Jul 2020 14:40:32 +0200 Subject: [PATCH 04/11] Update Utils.R --- R/Utils.R | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/R/Utils.R b/R/Utils.R index 2b5569b..8d60496 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -1712,7 +1712,11 @@ } anom = array(data = NA, dim = dim(data)) - clim = array(data = NA, dim = c(dim(data)['fyear'],dim(data)[member_dim])) + 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'])) + } for (i in 1:n_fyears){ if (member_dim %in% names(dim(data))){ for (m in 1:as.numeric(dim(data)[member_dim])){ -- GitLab From 142593d0e7f0d5727b0f10a865347932dafd727b Mon Sep 17 00:00:00 2001 From: cdelgado Date: Tue, 28 Jul 2020 16:12:17 +0200 Subject: [PATCH 05/11] Update AMV.R --- R/AMV.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/AMV.R b/R/AMV.R index 223edd2..90554f7 100644 --- a/R/AMV.R +++ b/R/AMV.R @@ -1,7 +1,7 @@ #'Compute the Atlantic Multidecadal Variability (AMV) index #' #'@description The Atlantic Multidecadal Variability (AMV), also known as Atlantic Multidecadal Oscillation (AMO), -#'is a model of natural variability of the sea surface temperatures (SST) over the North Atlantic Ocean +#'is a mode of natural variability of the sea surface temperatures (SST) over the North Atlantic Ocean #'on multi-decadal time scales. The AMV index is computed as the difference of weighted-averaged SST anomalies over #'the North Atlantic region (0ºN-60ºN, 280ºE-360ºE) and the weighted-averaged SST anomalies over 60ºS-60ºN, 0ºE-360ºE #'(Trenberth & Dennis, 2005; Doblas-Reyes et al., 2013). -- GitLab From 1d39125659697876e71fd320b5c46d16ee9e2f3e Mon Sep 17 00:00:00 2001 From: cdelgado Date: Mon, 24 Aug 2020 11:42:04 +0200 Subject: [PATCH 06/11] Update Utils.R --- R/Utils.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/Utils.R b/R/Utils.R index 8d60496..531a004 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -1751,7 +1751,7 @@ } else { target_dims = year_dim } - clim = multiApply::Apply(data = startR::Subset(x = data, along = year_dim, indices = indices_for_clim), + 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, -- GitLab From 7714e3f2ac5ec8f93d5343e14b552452aa88cd8c Mon Sep 17 00:00:00 2001 From: aho Date: Tue, 3 Nov 2020 12:01:21 +0100 Subject: [PATCH 07/11] Revise AMV() header, create unit tests. --- NAMESPACE | 9 ++ R/AMV.R | 273 ++++++++++++++++++++++++++------------ R/Utils.R | 16 ++- man/AMV.Rd | 117 ++++++++++++++++ man/GMST.Rd | 99 ++++++++++++++ man/GSAT.Rd | 81 +++++++++++ man/SPOD.Rd | 82 ++++++++++++ man/TPI.Rd | 81 +++++++++++ tests/testthat/test-AMV.R | 240 +++++++++++++++++++++++++++++++++ 9 files changed, 909 insertions(+), 89 deletions(-) create mode 100644 man/AMV.Rd create mode 100644 man/GMST.Rd create mode 100644 man/GSAT.Rd create mode 100644 man/SPOD.Rd create mode 100644 man/TPI.Rd create mode 100644 tests/testthat/test-AMV.R diff --git a/NAMESPACE b/NAMESPACE index 627bf51..b722dd3 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,6 @@ # Generated by roxygen2: do not edit by hand +export(AMV) export(AnimateMap) export(Clim) export(ColorBar) @@ -17,6 +18,8 @@ export(ConfigShowSimilarEntries) export(ConfigShowTable) export(Corr) export(Eno) +export(GMST) +export(GSAT) export(InsertDim) export(LeapYear) export(Load) @@ -31,11 +34,14 @@ export(RMS) export(RMSSS) export(Regression) export(Reorder) +export(SPOD) export(Season) +export(TPI) export(ToyModel) export(Trend) export(clim.colors) export(clim.palette) +import(ClimProjDiags) import(GEOmap) import(bigmemory) import(geomapdata) @@ -46,7 +52,10 @@ import(methods) import(multiApply) import(ncdf4) import(parallel) +import(s2dv) +import(startR) importFrom(ClimProjDiags,Subset) +importFrom(ClimProjDiags,WeightedMean) importFrom(abind,abind) importFrom(abind,adrop) importFrom(grDevices,bmp) diff --git a/R/AMV.R b/R/AMV.R index 90554f7..f500093 100644 --- a/R/AMV.R +++ b/R/AMV.R @@ -1,113 +1,217 @@ #'Compute the Atlantic Multidecadal Variability (AMV) index #' -#'@description The Atlantic Multidecadal Variability (AMV), also known as Atlantic Multidecadal Oscillation (AMO), -#'is a mode of natural variability of the sea surface temperatures (SST) over the North Atlantic Ocean -#'on multi-decadal time scales. The AMV index is computed as the difference of weighted-averaged SST anomalies over -#'the North Atlantic region (0ºN-60ºN, 280ºE-360ºE) and the weighted-averaged SST anomalies over 60ºS-60ºN, 0ºE-360ºE -#'(Trenberth & Dennis, 2005; Doblas-Reyes et al., 2013). -#'@author Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} -#'@author Roberto Bilbao, \email{roberto.bilbao@bsc.es} -#'@author Núria Pérez-Zanón, \email{nuria.perez@bsc.es} +#'The Atlantic Multidecadal Variability (AMV), also known as Atlantic +#'Multidecadal Oscillation (AMO), is a mode of natural variability of the sea +#'surface temperatures (SST) over the North Atlantic Ocean on multi-decadal +#'time scales. The AMV index is computed as the difference of weighted-averaged +#'SST anomalies over the North Atlantic region (0ºN-60ºN, 280ºE-360ºE) and the +#'weighted-averaged SST anomalies over 60ºS-60ºN, 0ºE-360ºE (Trenberth & +#'Dennis, 2005; Doblas-Reyes et al., 2013). #' -#'@param data Data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), -#' with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). -#' This data has to be provided, at least, over the whole region needed to compute the index. -#'@param data_lats An array with the latitudes of the data. -#'@param data_lons An array with the longitudes of the data. -#'@param type A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses). -#'@param lat_dim A string with the name of the latitude dimension ('lat' by default). -#'@param lon_dim A string with the name of the longitude dimension ('lon' by default). -#'@param mask An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a -#' combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the -#' observational dataset for a fair comparison between the forecast system and the reference dataset. -#'@param monini Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'. -#'@param fmonth_dim A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'. -#'@param sdate_dim A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'. -#'@param indices_for_clim Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). -#' In case of type='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. -#'@param year_dim A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'. -#'@param month_dim A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'. -#'@param member_dim A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'. +#'@param data A numerical array indicating the data to be used for the index +#' computation with latitude, longitude, start date, forecast month, and +#' member dimensions (in case of decadal predictions), with latitude, +#' longitude, year, month and member dimensions (in case of historical +#' simulations), or with latitude, longitude, year and month (in case of +#' observations or reanalyses). 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. +#'@param type A character string indicating the type of data ('dcpp' for +#' decadal predictions, 'hist' for historical simulations, or 'obs' for +#' observations or reanalyses). +#'@param lat_dim A character string of the name of the latitude dimension. The +#' default value is 'lat'. +#'@param lon_dim A character string of the name of the longitude dimension. The +#' default value is 'lon'. +#'@param mask An array of a mask (with 0's in the grid points that have to be +#' masked) or NULL (i.e., no mask is used). This parameter allows to remove +#' the values over land in case the dataset is a combination of surface air +#' temperature over land and sea surface temperature over the ocean. Also, it +#' can be used to mask those grid points that are missing in the observational +#' dataset for a fair comparison between the forecast system and the reference +#' dataset. The default value is NULL. +#'@param monini An integer indicating the month in which the forecast system is +#' initialized. Only used when parameter 'type' is 'dcpp'. The default value +#' is 11, i.e., initialized in November. +#'@param fmonth_dim A character string indicating the name of the forecast +#' month dimension. Only used if parameter 'type' is 'dcpp'. The default value +#' is 'fmonth'. +#'@param sdate_dim A character string indicating the name of the start date +#' dimension. Only used if parameter 'type' is 'dcpp'. The default value is +#' 'sdate'. +#'@param indices_for_clim A numeric vector of the indices of the years to +#' compute the climatology for calculating the anomalies, or NULL so the +#' climatology is calculated over the whole period. If the data are already +#' 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. +#'@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'. +#'@param month_dim A character string indicating the name of the month +#' dimension. The default value is 'month'. Only used if parameter 'type' is +#' 'hist' or 'obs'. +#'@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'. #' -#'@return The AMV index as funcion of the sdate, forecast year, and member (in case of decadal predictions); -#' as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +#'@return An numerical array of the AMV index with the dimensions of: +#' 1) sdate, forecast year, and member (in case of decadal predictions); +#' 2) year and the member (in case of historical simulations); or +#' 3) year (in case of observations or reanalyses). #' #'@examples #' ## Observations or reanalyses -#' obs = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) -#' lat = seq(-90, 90, 10) -#' lon = seq(0, 360, 10) -#' index_obs = AMV(data = obs, data_lats = lat, data_lons = lon, type = 'obs') +#' obs <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +#' lat <- seq(-90, 90, 10) +#' lon <- seq(0, 360, 10) +#' index_obs <- AMV(data = obs, data_lats = lat, data_lons = lon, type = 'obs') #' #' ## Historical simulations -#' hist = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) -#' lat = seq(-90, 90, 10) -#' lon = seq(0, 360, 10) -#' index_hist = AMV(data = hist, data_lats = lat, data_lons = lon, type = 'hist') +#' hist <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +#' lat <- seq(-90, 90, 10) +#' lon <- seq(0, 360, 10) +#' index_hist <- AMV(data = hist, data_lats = lat, data_lons = lon, type = 'hist') #' #' ## Decadal predictions -#' dcpp = array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) -#' lat = seq(-90, 90, 10) -#' lon = seq(0, 360, 10) -#' index_dcpp = AMV(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) +#' dcpp <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +#' lat <- seq(-90, 90, 10) +#' lon <- seq(0, 360, 10) +#' index_dcpp <- AMV(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) +#' +#'@author Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} +#'@author Roberto Bilbao, \email{roberto.bilbao@bsc.es} +#'@author Núria Pérez-Zanón, \email{nuria.perez@bsc.es} #' -#'@import ClimProjDiags +#'@importFrom ClimProjDiags WeightedMean #'@import multiApply -#'@import s2dv -#'@import startR #'@export 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'){ + indices_for_clim = NULL, year_dim = 'year', month_dim = 'month', + member_dim = 'member') { - ## Checkings - if (!is.array(data)){ - stop('data must be an array') + ## Input Checks + # data + if (is.null(data)) { + stop("Parameter 'data' cannot be NULL.") } - if(!class(data_lats)=='numeric'){ - stop('data_lats must be a numeric vector') + if (!is.numeric(data)) { + stop("Parameter 'data' must be a numeric array.") } - if (!class(data_lons)=='numeric'){ - stop('data_lons must be a numeric vector or NULL') + # data_lats and data_lons part1 + if(!(class(data_lats) == 'numeric' | class(data_lats) == 'integer')) { + stop("Parameter 'data_lats' must be a numeric vector.") } - if (!type %in% c('dcpp','hist','obs')){ - stop("type must be 'dcpp', 'hist' or 'obs'") + if(!(class(data_lons) == 'numeric' | class(data_lons) == 'integer')) { + stop("Parameter 'data_lons' must be a numeric vector.") } - if (!is.character(lat_dim)){ - stop('lat_dim must be a string') + # type + if (!type %in% c('dcpp', 'hist', 'obs')) { + stop("Parameter 'type' must be 'dcpp', 'hist', or 'obs'.") } - if (!is.character(lon_dim)){ - stop('lon_dim must be a string') + # lat_dim + if (!(is.character(lat_dim) & length(lat_dim) == 1)) { + stop("Parameter 'lat_dim' must be a character string.") } - if (!monini){ - stop("monini must be an integer from 1 to 12") + if (!lat_dim %in% names(dim(data))) { + stop("Parameter 'lat_dim' is not found in 'data' dimension.") } - if (!is.character(fmonth_dim)){ - stop('fmonth_dim must be a string') + # lon_dim + if (!(is.character(lon_dim) & length(lon_dim) == 1)) { + stop("Parameter 'lon_dim' must be a character string.") } - if (!is.character(sdate_dim)){ - stop('sdate_dim must be a string') + if (!lon_dim %in% names(dim(data))) { + stop("Parameter 'lon_dim' is not found in 'data' dimension.") } - if (!is.null(indices_for_clim) & !class(indices_for_clim) %in% c('numeric','integer') & !is.null(indices_for_clim) & !isFALSE(indices_for_clim)){ - stop("indices_for_clim must be a numeric vector, NULL to compute the anomalies based on the whole period, or FALSE if data are already anomalies") + # data_lats and data_lons part2 + if (dim(data)[lat_dim] != length(data_lats)){ + stop(paste0("The latitude dimension of parameter 'data' must be the same", + " length of parameter 'data_lats'.")) } - if (!is.character(year_dim)){ - stop('year_dim must be a string') + if (dim(data)[lon_dim] != length(data_lons)){ + stop(paste0("The longitude dimension of parameter 'data' must be the same", + " length of parameter 'data_lons'.")) } - if (!is.character(month_dim)){ - stop('month_dim must be a string') - } - if (!is.null(mask)){ - if (is.array(mask) & identical(names(dim(mask)),c(lat_dim,lon_dim)) & identical(as.integer(dim(mask)),c(length(data_lats),length(data_lons)))){ + # mask + if (!is.null(mask)) { + if (is.array(mask) & identical(names(dim(mask)), c(lat_dim,lon_dim)) & + identical(as.integer(dim(mask)), c(length(data_lats), length(data_lons)))) { ## To mask those grid point that are missing in the observations - mask <- s2dv::Reorder(data = mask, order = c(lat_dim,lon_dim)) + mask <- s2dv::Reorder(data = mask, order = c(lat_dim, lon_dim)) fun_mask <- function(data, mask){ data[mask == 0] <- NA return(data) } - data <- multiApply::Apply(data = data, target_dims = c(lat_dim,lon_dim), fun = fun_mask, mask = mask)$output1 + data <- multiApply::Apply(data = data, target_dims = c(lat_dim, lon_dim), + fun = fun_mask, mask = mask)$output1 } else { - stop('mask must be NULL (no mask) or an array with c(lat_dim,lon_dim) dimensions and 0 in those grid points that have to be masked') + stop(paste0("Parameter 'mask' must be NULL (no mask) or a numerical array ", + "with c(lat_dim, lon_dim) dimensions and 0 in those grid ", + "points that have to be masked.")) + } + } + # monini + if (type == 'dcpp') { + if (!is.numeric(monini) | monini %% 1 != 0 | monini < 1 | + length(monini) > 12) { + stop("Parameter 'monini' must be an integer from 1 to 12.") + } + } + # fmonth_dim + if (type == 'dcpp') { + if (!(is.character(fmonth_dim) & length(fmonth_dim) == 1)) { + stop("Parameter 'fmonth_dim' must be a character string.") + } + if (!fmonth_dim %in% names(dim(data))) { + stop("Parameter 'fmonth_dim' is not found in 'data' dimension.") + } + } + # sdate_dim + if (type == 'dcpp') { + if (!(is.character(sdate_dim) & length(sdate_dim) == 1)) { + stop("Parameter 'sdate_dim' must be a character string.") + } + if (!sdate_dim %in% names(dim(data))) { + stop("Parameter 'sdate_dim' is not found in 'data' dimension.") + } + } + # indices_for_clim + if (!is.null(indices_for_clim)) { + if (!class(indices_for_clim) %in% c('numeric', 'integer') + & !(is.logical(indices_for_clim) & !any(indices_for_clim))) { + stop(paste0("The parameter 'indices_for_clim' must be a numeric vector ", + "or NULL to compute the anomalies based on the whole period, ", + "or FALSE if data are already anomalies")) + } + } + # year_dim + if (type == 'hist' | type == 'obs') { + if (!(is.character(year_dim) & length(year_dim) == 1)) { + stop("Parameter 'year_dim' must be a character string.") + } + if (!year_dim %in% names(dim(data))) { + stop("Parameter 'year_dim' is not found in 'data' dimension.") + } + } + # month_dim + if (type == 'hist' | type == 'obs') { + if (!(is.character(month_dim) & length(month_dim) == 1)) { + stop("Parameter 'month_dim' must be a character string.") + } + if (!month_dim %in% names(dim(data))) { + 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.") } } @@ -120,15 +224,20 @@ AMV <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'lo regions$reg1 <- c(lon_min_1, lon_max_1, lat_min_1, lat_max_1) regions$reg2 <- c(lon_min_2, lon_max_2, lat_min_2, lat_max_2) - mean_1 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = regions$reg1, - londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) - mean_2 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = regions$reg2, - londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) - data <- ClimProjDiags::CombineIndices(indices = list(mean_1,mean_2), weights = NULL, operation = 'subtract') # (mean_1 - mean_2) + mean_1 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, + region = regions$reg1, + londim = which(names(dim(data)) == lon_dim), + latdim = which(names(dim(data)) == lat_dim)) + mean_2 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, + region = regions$reg2, + londim = which(names(dim(data)) == lon_dim), + latdim = which(names(dim(data)) == lat_dim)) + 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) return(INDEX) -} \ No newline at end of file +} diff --git a/R/Utils.R b/R/Utils.R index 531a004..dd6b898 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -1686,10 +1686,11 @@ data = s2dv::Reorder(data = data, order = c('fyear',sdate_dim)) } - if (isFALSE(indices_for_clim)){ - + if (is.logical(indices_for_clim)) { + if(!any(indices_for_clim)) { # indices_for_clim == FALSE -> anomalies are directly given - anom = data + anom = data + } } else { @@ -1734,10 +1735,11 @@ data = multiApply::Apply(data = data, target_dims = month_dim, fun = mean)$output1 - if (isFALSE(indices_for_clim)){ - - anom = data - + if (is.logical(indices_for_clim)) { + if(!any(indices_for_clim)) { + anom = data + } + } else { if (is.null(indices_for_clim)){ diff --git a/man/AMV.Rd b/man/AMV.Rd new file mode 100644 index 0000000..3e87dc7 --- /dev/null +++ b/man/AMV.Rd @@ -0,0 +1,117 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AMV.R +\name{AMV} +\alias{AMV} +\title{Compute the Atlantic Multidecadal Variability (AMV) index} +\usage{ +AMV(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") +} +\arguments{ +\item{data}{A numerical array indicating the data to be used for the index +computation with latitude, longitude, start date, forecast month, and +member dimensions (in case of decadal predictions), with latitude, +longitude, year, month and member dimensions (in case of historical +simulations), or with latitude, longitude, year and month (in case of +observations or reanalyses). This data has to be provided, at least, over +the whole region needed to compute the index.} + +\item{data_lats}{A numeric vector indicating the latitudes of the data.} + +\item{data_lons}{A numeric vector indicating the longitudes of the data.} + +\item{type}{A character string indicating the type of data ('dcpp' for +decadal predictions, 'hist' for historical simulations, or 'obs' for +observations or reanalyses).} + +\item{lat_dim}{A character string of the name of the latitude dimension. The +default value is 'lat'.} + +\item{lon_dim}{A character string of the name of the longitude dimension. The +default value is 'lon'.} + +\item{mask}{An array of a mask (with 0's in the grid points that have to be +masked) or NULL (i.e., no mask is used). This parameter allows to remove +the values over land in case the dataset is a combination of surface air +temperature over land and sea surface temperature over the ocean. Also, it +can be used to mask those grid points that are missing in the observational +dataset for a fair comparison between the forecast system and the reference +dataset. The default value is NULL.} + +\item{monini}{An integer indicating the month in which the forecast system is +initialized. Only used when parameter 'type' is 'dcpp'. The default value +is 11, i.e., initialized in November.} + +\item{fmonth_dim}{A character string indicating the name of the forecast +month dimension. Only used if parameter 'type' is 'dcpp'. The default value +is 'fmonth'.} + +\item{sdate_dim}{A character string indicating the name of the start date +dimension. Only used if parameter 'type' is 'dcpp'. The default value is +'sdate'.} + +\item{indices_for_clim}{A numeric vector of the indices of the years to +compute the climatology for calculating the anomalies, or NULL so the +climatology is calculated over the whole period. If the data are already +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.} + +\item{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'.} + +\item{month_dim}{A character string indicating the name of the month +dimension. The default value is 'month'. Only used if parameter 'type' is +'hist' or 'obs'.} + +\item{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'.} +} +\value{ +An numerical array of the AMV index with the dimensions of: + 1) sdate, forecast year, and member (in case of decadal predictions); + 2) year and the member (in case of historical simulations); or + 3) year (in case of observations or reanalyses). +} +\description{ +The Atlantic Multidecadal Variability (AMV), also known as Atlantic +Multidecadal Oscillation (AMO), is a mode of natural variability of the sea +surface temperatures (SST) over the North Atlantic Ocean on multi-decadal +time scales. The AMV index is computed as the difference of weighted-averaged +SST anomalies over the North Atlantic region (0ºN-60ºN, 280ºE-360ºE) and the +weighted-averaged SST anomalies over 60ºS-60ºN, 0ºE-360ºE (Trenberth & +Dennis, 2005; Doblas-Reyes et al., 2013). +} +\examples{ +## Observations or reanalyses +obs <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +lat <- seq(-90, 90, 10) +lon <- seq(0, 360, 10) +index_obs <- AMV(data = obs, data_lats = lat, data_lons = lon, type = 'obs') + +## Historical simulations +hist <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +lat <- seq(-90, 90, 10) +lon <- seq(0, 360, 10) +index_hist <- AMV(data = hist, data_lats = lat, data_lons = lon, type = 'hist') + +## Decadal predictions +dcpp <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +lat <- seq(-90, 90, 10) +lon <- seq(0, 360, 10) +index_dcpp <- AMV(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) + +} +\author{ +Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} + +Roberto Bilbao, \email{roberto.bilbao@bsc.es} + +Núria Pérez-Zanón, \email{nuria.perez@bsc.es} +} + diff --git a/man/GMST.Rd b/man/GMST.Rd new file mode 100644 index 0000000..78508b5 --- /dev/null +++ b/man/GMST.Rd @@ -0,0 +1,99 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/GMST.R +\name{GMST} +\alias{GMST} +\title{Compute the Global Mean Surface Temperature (GMST) anomalies} +\usage{ +GMST(data_tas, data_tos, data_lats, data_lons, mask_sea_land, sea_value, type, + mask = NULL, lat_dim = "lat", lon_dim = "lon", monini = 11, + fmonth_dim = "fmonth", sdate_dim = "sdate", indices_for_clim = NULL, + year_dim = "year", month_dim = "month", member_dim = "member") +} +\arguments{ +\item{data_tas}{Surface air temperature data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), +with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). +This data has to be provided, at least, over the whole region needed to compute the index. The dimensions must be identical to those of data_tos.} + +\item{data_tos}{Sea surface temperature data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), +with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). +This data has to be provided, at least, over the whole region needed to compute the index. The dimensions must be identical to those of data_tas.} + +\item{data_lats}{An array with the latitudes of the data.} + +\item{data_lons}{An array with the longitudes of the data.} + +\item{mask_sea_land}{# mask for blending the surface air temperature over land and the sea surface temperature over the ocean. It must have lat_dim and lon_dim dimensions, and their lengths +must be equal to the length of data_lats and data_lons.} + +\item{sea_value}{Value of the sea grid points in the mask.} + +\item{type}{A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses).} + +\item{mask}{An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a +combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the +observational dataset for a fair comparison between the forecast system and the reference dataset.} + +\item{lat_dim}{A string with the name of the latitude dimension ('lat' by default).} + +\item{lon_dim}{A string with the name of the longitude dimension ('lon' by default).} + +\item{monini}{Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'.} + +\item{fmonth_dim}{A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'.} + +\item{sdate_dim}{A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'.} + +\item{indices_for_clim}{Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). +In case of type='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.} + +\item{year_dim}{A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'.} + +\item{month_dim}{A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'.} + +\item{member_dim}{A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'.} +} +\value{ +The GMST anomalies as function of the sdate, forecast year, and member (in case of decadal predictions); + as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +} +\description{ +The Global Mean Surface Temperature (GMST) anomalies are computed as the weighted-averaged +surface air temperature anomalies over land and sea surface temperature anomalies over the ocean. +} +\examples{ +## Observations or reanalyses +obs_tas = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +obs_tos = array(2:101, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +mask_sea_land = array(c(1,0,1), dim = c(lat = 19, lon = 37)) +sea_value = 1 +lat = seq(-90, 90, 10) +lon = seq(0, 360, 10) +index_obs = GMST(data_tas = obs_tas, data_tos = obs_tos, data_lats = lat, data_lons = lon, type = 'obs', mask_sea_land = mask_sea_land, sea_value = sea_value) + +## Historical simulations +hist_tas = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +hist_tos = array(2:101, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +mask_sea_land = array(c(1,0,1), dim = c(lat = 19, lon = 37)) +sea_value = 1 +lat = seq(-90, 90, 10) +lon = seq(0, 360, 10) +index_hist = GMST(data_tas = hist_tas, data_tos = hist_tos, data_lats = lat, data_lons = lon, type = 'hist', mask_sea_land = mask_sea_land, sea_value = sea_value) + +## Decadal predictions +dcpp_tas = array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +dcpp_tos = array(2:101, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +mask_sea_land = array(c(1,0,1), dim = c(lat = 19, lon = 37)) +sea_value = 1 +lat = seq(-90, 90, 10) +lon = seq(0, 360, 10) +index_dcpp = GMST(data_tas = dcpp_tas, data_tos = dcpp_tos, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1, mask_sea_land = mask_sea_land, sea_value = sea_value) + +} +\author{ +Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} + +Roberto Bilbao, \email{roberto.bilbao@bsc.es} + +Núria Pérez-Zanón, \email{nuria.perez@bsc.es} +} + diff --git a/man/GSAT.Rd b/man/GSAT.Rd new file mode 100644 index 0000000..fda763e --- /dev/null +++ b/man/GSAT.Rd @@ -0,0 +1,81 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/GSAT.R +\name{GSAT} +\alias{GSAT} +\title{Compute the Global Surface Air Temperature (GSAT) anomalies} +\usage{ +GSAT(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") +} +\arguments{ +\item{data}{Data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), +with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). +This data has to be provided, at least, over the whole region needed to compute the index.} + +\item{data_lats}{An array with the latitudes of the data.} + +\item{data_lons}{An array with the longitudes of the data.} + +\item{type}{A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses).} + +\item{lat_dim}{A string with the name of the latitude dimension ('lat' by default).} + +\item{lon_dim}{A string with the name of the longitude dimension ('lon' by default).} + +\item{mask}{An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a +combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the +observational dataset for a fair comparison between the forecast system and the reference dataset.} + +\item{monini}{Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'.} + +\item{fmonth_dim}{A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'.} + +\item{sdate_dim}{A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'.} + +\item{indices_for_clim}{Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). +In case of type='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.} + +\item{year_dim}{A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'.} + +\item{month_dim}{A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'.} + +\item{member_dim}{A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'.} +} +\value{ +The GSAT anomalies as function of the sdate, forecast year, and member (in case of decadal predictions); + as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +} +\description{ +The Global Surface Air Temperature (GSAT) anomalies are computed as the +weighted-averaged surface air temperature anomalies over the global region. +} +\examples{ +## Observations or reanalyses +obs = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +lat = seq(-90, 90, 10) +lon = seq(0, 360, 10) +index_obs = GSAT(data = obs, data_lats = lat, data_lons = lon, type = 'obs') + +## Historical simulations +hist = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +lat = seq(-90, 90, 10) +lon = seq(0, 360, 10) +index_hist = GSAT(data = hist, data_lats = lat, data_lons = lon, type = 'hist') + +## Decadal predictions +dcpp = array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +lat = seq(-90, 90, 10) +lon = seq(0, 360, 10) +index_dcpp = GSAT(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) + +} +\author{ +Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} + +Núria Pérez-Zanón, \email{nuria.perez@bsc.es} + +Roberto Bilbao, \email{roberto.bilbao@bsc.es} +} + diff --git a/man/SPOD.Rd b/man/SPOD.Rd new file mode 100644 index 0000000..36ad341 --- /dev/null +++ b/man/SPOD.Rd @@ -0,0 +1,82 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/SPOD.R +\name{SPOD} +\alias{SPOD} +\title{Compute the South Pacific Ocean Dipole (SPOD) index} +\usage{ +SPOD(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") +} +\arguments{ +\item{data}{Data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), +with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). +This data has to be provided, at least, over the whole region needed to compute the index.} + +\item{data_lats}{An array with the latitudes of the data.} + +\item{data_lons}{An array with the longitudes of the data.} + +\item{type}{A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses).} + +\item{lat_dim}{A string with the name of the latitude dimension ('lat' by default).} + +\item{lon_dim}{A string with the name of the longitude dimension ('lon' by default).} + +\item{mask}{An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a +combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the +observational dataset for a fair comparison between the forecast system and the reference dataset.} + +\item{monini}{Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'.} + +\item{fmonth_dim}{A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'.} + +\item{sdate_dim}{A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'.} + +\item{indices_for_clim}{Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). +In case of type='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.} + +\item{year_dim}{A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'.} + +\item{month_dim}{A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'.} + +\item{member_dim}{A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'.} +} +\value{ +The SPOD index as function of the sdate, forecast year, and member (in case of decadal predictions); + as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +} +\description{ +The South Pacific Ocean Dipole (SPOD) index is related to the El Nino-Southern Oscillation (ENSO) and the Inderdecadal Pacific Oscillation (IPO). +The SPOD index is computed as the difference of weighted-averaged SST anomalies over 20ºS-48ºS, 165ºE-190ºE (NW pole) and the weighted-averaged SST anomalies over +44ºS-65ºS, 220ºE-260ºE (SE pole) (Saurral et al., 2020). +} +\examples{ +## Observations or reanalyses +obs = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +lat = seq(-90, 90, 10) +lon = seq(0, 360, 10) +index_obs = SPOD(data = obs, data_lats = lat, data_lons = lon, type = 'obs') + +## Historical simulations +hist = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +lat = seq(-90, 90, 10) +lon = seq(0, 360, 10) +index_hist = SPOD(data = hist, data_lats = lat, data_lons = lon, type = 'hist') + +## Decadal predictions +dcpp = array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +lat = seq(-90, 90, 10) +lon = seq(0, 360, 10) +index_dcpp = SPOD(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) + +} +\author{ +Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} + +Núria Pérez-Zanón, \email{nuria.perez@bsc.es} + +Roberto Bilbao, \email{roberto.bilbao@bsc.es} +} + diff --git a/man/TPI.Rd b/man/TPI.Rd new file mode 100644 index 0000000..0e5655b --- /dev/null +++ b/man/TPI.Rd @@ -0,0 +1,81 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/TPI.R +\name{TPI} +\alias{TPI} +\title{Compute the Tripole Index (TPI) for the Interdecadal Pacific Oscillation (IPO)} +\usage{ +TPI(data, data_lats, data_lons, index, 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") +} +\arguments{ +\item{data}{Data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), +with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). +This data has to be provided, at least, over the whole region needed to compute the index.} + +\item{data_lats}{An array with the latitudes of the data.} + +\item{data_lons}{An array with the longitudes of the data.} + +\item{type}{A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses).} + +\item{lat_dim}{A string with the name of the latitude dimension ('lat' by default).} + +\item{lon_dim}{A string with the name of the longitude dimension ('lon' by default).} + +\item{mask}{An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a +combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the +observational dataset for a fair comparison between the forecast system and the reference dataset.} + +\item{monini}{Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'.} + +\item{fmonth_dim}{A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'.} + +\item{sdate_dim}{A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'.} + +\item{indices_for_clim}{Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). +In case of type='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.} + +\item{year_dim}{A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'.} + +\item{month_dim}{A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'.} + +\item{member_dim}{A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'.} +} +\value{ +The TPI as function of the sdate, forecast year, and member (in case of decadal predictions); + as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +} +\description{ +The Tripole Index (TPI) for the Interdecadal Pacific Oscillation (IPO) is computed as the difference of weighted-averaged SST anomalies over +10ºS-10ºN, 170ºE-270ºE minus the mean of the weighted-averaged SST anomalies over 25ºN-45ºN, 140ºE-215ºE and 50ºS-15ºS, 150ºE-200ºE (Henley et al., 2015). +} +\examples{ +## Observations or reanalyses +obs = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +lat = seq(-90, 90, 10) +lon = seq(0, 360, 10) +index_obs = TPI(data = obs, data_lats = lat, data_lons = lon, type = 'obs') + +## Historical simulations +hist = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +lat = seq(-90, 90, 10) +lon = seq(0, 360, 10) +index_hist = TPI(data = hist, data_lats = lat, data_lons = lon, type = 'hist') + +## Decadal predictions +dcpp = array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +lat = seq(-90, 90, 10) +lon = seq(0, 360, 10) +index_dcpp = TPI(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) + +} +\author{ +Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} + +Núria Pérez-Zanón, \email{nuria.perez@bsc.es} + +Roberto Bilbao, \email{roberto.bilbao@bsc.es} +} + diff --git a/tests/testthat/test-AMV.R b/tests/testthat/test-AMV.R new file mode 100644 index 0000000..562605b --- /dev/null +++ b/tests/testthat/test-AMV.R @@ -0,0 +1,240 @@ +context("s2dv::AMV tests") + +############################################## + lat <- seq(-90, 90, 10) + lon <- seq(0, 360, 10) + set.seed(1) + mask <- array(sample(c(0, 1), 19*37, replace = T), dim = c(lat = 19, lon = 37)) + +## Observations or reanalyses + dat1 <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) + + ## Historical simulations + dat2 <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) + + ## Decadal predictions + dat3 <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) + +############################################## +test_that("1. Input checks", { + + expect_error( + AMV(data = c()), + "Parameter 'data' cannot be NULL." + ) + expect_error( + AMV(data = 'a'), + "Parameter 'data' must be a numeric array." + ) + expect_error( + AMV(data = dat1, data_lats = c()), + "Parameter 'data_lats' must be a numeric vector." + ) + expect_error( + AMV(data = dat1, data_lats = lat, data_lons = c()), + "Parameter 'data_lons' must be a numeric vector." + ) + expect_error( + AMV(data = dat1, data_lats = lat, data_lons = lon, type = TRUE), + "Parameter 'type' must be 'dcpp', 'hist', or 'obs'." + ) + expect_error( + AMV(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + lat_dim = 2), + "Parameter 'lat_dim' must be a character string." + ) + expect_error( + AMV(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + lat_dim = 'lats'), + "Parameter 'lat_dim' is not found in 'data' dimension." + ) + expect_error( + AMV(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + lon_dim = 2), + "Parameter 'lon_dim' must be a character string." + ) + expect_error( + AMV(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + lon_dim = 'lons'), + "Parameter 'lon_dim' is not found in 'data' dimension." + ) + expect_error( + AMV(data = dat3, data_lats = seq(-90, 90, 11), data_lons = lon, type = 'dcpp'), + paste0("The latitude dimension of parameter 'data' must be the same", + " length of parameter 'data_lats'.") + ) + expect_error( + AMV(data = dat3, data_lats = lat, data_lons = c(1:10), type = 'dcpp'), + paste0("The longitude dimension of parameter 'data' must be the same", + " length of parameter 'data_lons'.") + ) + expect_error( + AMV(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + mask = c(1, 0, 1, 0, 0)), + paste0("Parameter 'mask' must be NULL \\(no mask\\) or a numerical array ", + "with c\\(lat_dim, lon_dim\\) dimensions and 0 in those grid ", + "points that have to be masked.") + ) + expect_error( + AMV(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + mask = array(rep(c(1, 0), 100), dim = c(lat = 20, lon = 5))), + paste0("Parameter 'mask' must be NULL \\(no mask\\) or a numerical array ", + "with c\\(lat_dim, lon_dim\\) dimensions and 0 in those grid ", + "points that have to be masked.") + ) + expect_error( + AMV(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + monini = 0), + "Parameter 'monini' must be an integer from 1 to 12." + ) + expect_error( + AMV(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + fmonth_dim = c('lat', 'lon')), + "Parameter 'fmonth_dim' must be a character string." + ) + expect_error( + AMV(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + fmonth_dim = 'asd'), + "Parameter 'fmonth_dim' is not found in 'data' dimension." + ) + expect_error( + AMV(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + sdate_dim = 2), + "Parameter 'sdate_dim' must be a character string." + ) + expect_error( + AMV(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + sdate_dim = 'date'), + "Parameter 'sdate_dim' is not found in 'data' dimension." + ) + expect_error( + AMV(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + indices_for_clim = TRUE), + paste0("The parameter 'indices_for_clim' must be a numeric vector ", + "or NULL to compute the anomalies based on the whole period, ", + "or FALSE if data are already anomalies") + ) + expect_error( + AMV(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + year_dim = 1), + "Parameter 'year_dim' must be a character string." + ) + expect_error( + AMV(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + year_dim = 'yr'), + "Parameter 'year_dim' is not found in 'data' dimension." + ) + expect_error( + AMV(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + month_dim = 0), + "Parameter 'month_dim' must be a character string." + ) + expect_error( + AMV(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + month_dim = 'mth'), + "Parameter 'month_dim' is not found in 'data' dimension." + ) + expect_error( + AMV(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + member_dim = 1), + "Parameter 'member_dim' must be a character string." + ) + expect_error( + AMV(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + member_dim = 'ens'), + "Parameter 'member_dim' is not found in 'data' dimension." + ) + +}) + +############################################## +test_that("2. Output checks: dat1", { + res1_1 <- AMV(data = dat1, data_lats = lat, data_lons = lon, type = 'obs') + + expect_equal( + dim(res1_1), + c(year = 5) + ) + expect_equal( + mean(res1_1)*10^17, + -4.440892, + tolerance = 0.00001 + ) + + res1_2 <- AMV(data = dat1, data_lats = lat, data_lons = lon, type = 'obs', + indices_for_clim = F) + expect_equal( + mean(res1_2), + -1.78994, + tolerance = 0.00001 + ) + +}) +############################################## +test_that("3. Output checks: dat2", { + res2_1 <- AMV(data = dat2, data_lats = lat, data_lons = lon, type = 'hist') + expect_equal( + dim(res2_1), + c(year = 5, member = 5) + ) + expect_equal( + mean(res2_1)*10^17, + -3.549383, + tolerance = 0.00001 + ) + + res2_2 <- AMV(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + indices_for_clim = 2:3) + expect_equal( + mean(res2_2)*10^16, + 8.616545, + tolerance = 0.00001 + ) + + res2_3 <- AMV(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + mask = mask) + expect_equal( + mean(res2_3)*10^17, + -4.437162, + tolerance = 0.00001 + ) + +}) +############################################## +test_that("4. Output checks: dat3", { + res3_1 <- AMV(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp') + expect_equal( + dim(res3_1), + c(fyear = 1, sdate = 5, member = 5) + ) + expect_equal( + mean(res3_1)*10^18, + -8.881784, + tolerance = 0.00001 + ) + + res3_2 <- AMV(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + indices_for_clim = 2:3) + expect_equal( + mean(res3_2)*10^16, + 8.459899, + tolerance = 0.00001 + ) + + res3_3 <- AMV(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + monini = 5) + expect_equal( + mean(res3_3)*10^17, + -4.884981, + tolerance = 0.00001 + ) + + res3_4 <- AMV(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + mask = mask) + expect_equal( + mean(res3_4)*10^17, + 2.220446, + tolerance = 0.00001 + ) + +}) -- GitLab From 38dbb6b35759e2a98b3ae17329d51c51f0f28614 Mon Sep 17 00:00:00 2001 From: aho Date: Tue, 3 Nov 2020 12:12:26 +0100 Subject: [PATCH 08/11] Remove unecessary import packages --- NAMESPACE | 2 -- R/GMST.R | 4 +--- R/GSAT.R | 4 +--- R/SPOD.R | 4 +--- R/TPI.R | 4 +--- 5 files changed, 4 insertions(+), 14 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index b722dd3..9d15ffe 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -52,8 +52,6 @@ import(methods) import(multiApply) import(ncdf4) import(parallel) -import(s2dv) -import(startR) importFrom(ClimProjDiags,Subset) importFrom(ClimProjDiags,WeightedMean) importFrom(abind,abind) diff --git a/R/GMST.R b/R/GMST.R index 6ea8216..9c147be 100644 --- a/R/GMST.R +++ b/R/GMST.R @@ -65,8 +65,6 @@ #' #'@import ClimProjDiags #'@import multiApply -#'@import s2dv -#'@import startR #'@export GMST <- function(data_tas, data_tos, data_lats, data_lons, mask_sea_land, sea_value, type, mask = NULL, lat_dim = 'lat', lon_dim = 'lon', monini = 11, @@ -157,4 +155,4 @@ GMST <- function(data_tas, data_tos, data_lats, data_lons, mask_sea_land, sea_va sdate_dim = sdate_dim, year_dim = year_dim, month_dim = month_dim, member_dim = member_dim) return(INDEX) -} \ No newline at end of file +} diff --git a/R/GSAT.R b/R/GSAT.R index 8f8d4ab..ced384c 100644 --- a/R/GSAT.R +++ b/R/GSAT.R @@ -50,8 +50,6 @@ #' #'@import ClimProjDiags #'@import multiApply -#'@import s2dv -#'@import startR #'@export GSAT <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'lon', mask = NULL, monini = 11, fmonth_dim = 'fmonth', sdate_dim = 'sdate', @@ -116,4 +114,4 @@ GSAT <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'l sdate_dim = sdate_dim, year_dim = year_dim, month_dim = month_dim, member_dim = member_dim) return(INDEX) -} \ No newline at end of file +} diff --git a/R/SPOD.R b/R/SPOD.R index afc5a73..688a527 100644 --- a/R/SPOD.R +++ b/R/SPOD.R @@ -51,8 +51,6 @@ #' #'@import ClimProjDiags #'@import multiApply -#'@import s2dv -#'@import startR #'@export SPOD <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'lon', mask = NULL, monini = 11, fmonth_dim = 'fmonth', sdate_dim = 'sdate', @@ -130,4 +128,4 @@ SPOD <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'l sdate_dim = sdate_dim, year_dim = year_dim, month_dim = month_dim, member_dim = member_dim) return(INDEX) -} \ No newline at end of file +} diff --git a/R/TPI.R b/R/TPI.R index 4e6537e..fbf6e59 100644 --- a/R/TPI.R +++ b/R/TPI.R @@ -50,8 +50,6 @@ #' #'@import ClimProjDiags #'@import multiApply -#'@import s2dv -#'@import startR #'@export TPI <- function(data, data_lats, data_lons, index, type, lat_dim = 'lat', lon_dim = 'lon', mask = NULL, monini = 11, fmonth_dim = 'fmonth', sdate_dim = 'sdate', @@ -136,4 +134,4 @@ TPI <- function(data, data_lats, data_lons, index, type, lat_dim = 'lat', lon_di sdate_dim = sdate_dim, year_dim = year_dim, month_dim = month_dim, member_dim = member_dim) return(INDEX) -} \ No newline at end of file +} -- GitLab From fd354f0494de859459da9c58285e711ff4191cb2 Mon Sep 17 00:00:00 2001 From: aho Date: Tue, 3 Nov 2020 14:01:16 +0100 Subject: [PATCH 09/11] Correct the monini check --- R/AMV.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/AMV.R b/R/AMV.R index f500093..5f8a824 100644 --- a/R/AMV.R +++ b/R/AMV.R @@ -156,7 +156,7 @@ AMV <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'lo # monini if (type == 'dcpp') { if (!is.numeric(monini) | monini %% 1 != 0 | monini < 1 | - length(monini) > 12) { + monini > 12) { stop("Parameter 'monini' must be an integer from 1 to 12.") } } -- GitLab From 5ca91e3f95bad42674771e01c406e8a7d901137c Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 4 Nov 2020 13:05:25 +0100 Subject: [PATCH 10/11] Review indices functions and create unit tests --- NAMESPACE | 2 +- R/AMV.R | 21 ++- R/GMST.R | 331 +++++++++++++++++++++++++------------ R/GSAT.R | 251 ++++++++++++++++++++-------- R/SPOD.R | 276 +++++++++++++++++++++---------- R/TPI.R | 249 ++++++++++++++++++++-------- man/AMV.Rd | 17 +- man/GMST.Rd | 171 +++++++++++-------- man/GSAT.Rd | 122 ++++++++------ man/SPOD.Rd | 130 +++++++++------ man/TPI.Rd | 104 ++++++++---- tests/testthat/test-GMST.R | 287 ++++++++++++++++++++++++++++++++ tests/testthat/test-GSAT.R | 240 +++++++++++++++++++++++++++ tests/testthat/test-SPOD.R | 240 +++++++++++++++++++++++++++ tests/testthat/test-TPI.R | 245 +++++++++++++++++++++++++++ 15 files changed, 2139 insertions(+), 547 deletions(-) create mode 100644 tests/testthat/test-GMST.R create mode 100644 tests/testthat/test-GSAT.R create mode 100644 tests/testthat/test-SPOD.R create mode 100644 tests/testthat/test-TPI.R diff --git a/NAMESPACE b/NAMESPACE index 9d15ffe..98c21bd 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -41,7 +41,6 @@ export(ToyModel) export(Trend) export(clim.colors) export(clim.palette) -import(ClimProjDiags) import(GEOmap) import(bigmemory) import(geomapdata) @@ -52,6 +51,7 @@ import(methods) import(multiApply) import(ncdf4) import(parallel) +importFrom(ClimProjDiags,CombineIndices) importFrom(ClimProjDiags,Subset) importFrom(ClimProjDiags,WeightedMean) importFrom(abind,abind) diff --git a/R/AMV.R b/R/AMV.R index 5f8a824..6657f6f 100644 --- a/R/AMV.R +++ b/R/AMV.R @@ -8,13 +8,12 @@ #'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 indicating the data to be used for the index -#' computation with latitude, longitude, start date, forecast month, and -#' member dimensions (in case of decadal predictions), with latitude, -#' longitude, year, month and member dimensions (in case of historical -#' simulations), or with latitude, longitude, year and month (in case of -#' observations or reanalyses). This data has to be provided, at least, over -#' the whole region needed to compute the index. +#'@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 +#' 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. #'@param type A character string indicating the type of data ('dcpp' for @@ -57,9 +56,9 @@ #' dimension. The default value is 'member'. Only used if parameter 'type' is #' 'dcpp' or 'hist'. #' -#'@return An numerical array of the AMV index with the dimensions of: +#'@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 the member (in case of historical simulations); or +#' 2) year and member (in case of historical simulations); or #' 3) year (in case of observations or reanalyses). #' #'@examples @@ -102,10 +101,10 @@ AMV <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'lo stop("Parameter 'data' must be a numeric array.") } # data_lats and data_lons part1 - if(!(class(data_lats) == 'numeric' | class(data_lats) == 'integer')) { + if (!(class(data_lats) == 'numeric' | class(data_lats) == 'integer')) { stop("Parameter 'data_lats' must be a numeric vector.") } - if(!(class(data_lons) == 'numeric' | class(data_lons) == 'integer')) { + if (!(class(data_lons) == 'numeric' | class(data_lons) == 'integer')) { stop("Parameter 'data_lons' must be a numeric vector.") } # type diff --git a/R/GMST.R b/R/GMST.R index 9c147be..c827951 100644 --- a/R/GMST.R +++ b/R/GMST.R @@ -1,154 +1,275 @@ #'Compute the Global Mean Surface Temperature (GMST) anomalies #' -#'@description The Global Mean Surface Temperature (GMST) anomalies are computed as the weighted-averaged -#'surface air temperature anomalies over land and sea surface temperature anomalies over the ocean. -#'@author Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} -#'@author Roberto Bilbao, \email{roberto.bilbao@bsc.es} -#'@author Núria Pérez-Zanón, \email{nuria.perez@bsc.es} +#'The Global Mean Surface Temperature (GMST) anomalies are computed as the +#'weighted-averaged surface air temperature anomalies over land and sea surface +#'temperature anomalies over the ocean. #' -#'@param data_tas Surface air temperature data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), -#' with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). -#' This data has to be provided, at least, over the whole region needed to compute the index. The dimensions must be identical to those of data_tos. -#'@param data_tos Sea surface temperature data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), -#' with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). -#' This data has to be provided, at least, over the whole region needed to compute the index. The dimensions must be identical to those of data_tas. -#'@param data_lats An array with the latitudes of the data. -#'@param data_lons An array with the longitudes of the data. -#'@param mask_sea_land # mask for blending the surface air temperature over land and the sea surface temperature over the ocean. It must have lat_dim and lon_dim dimensions, and their lengths -#' must be equal to the length of data_lats and data_lons. -#'@param sea_value Value of the sea grid points in the mask. -#'@param type A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses). -#'@param lat_dim A string with the name of the latitude dimension ('lat' by default). -#'@param lon_dim A string with the name of the longitude dimension ('lon' by default). -#'@param mask An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a -#' combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the -#' observational dataset for a fair comparison between the forecast system and the reference dataset. -#'@param monini Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'. -#'@param fmonth_dim A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'. -#'@param sdate_dim A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'. -#'@param indices_for_clim Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). -#' In case of type='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. -#'@param year_dim A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'. -#'@param month_dim A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'. -#'@param member_dim A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'. +#'@param data_tas A numerical array indicating the surface air temperature data +#' 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 provided, at least, +#' over the whole region needed to compute the index. The dimensions must be +#' identical to those of data_tos. +#'@param data_tos A numerical array indicating the sea surface temperature data +#' 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 provided, at least, +#' over the whole region needed to compute the index. The dimensions must be +#' identical to those of data_tas. +#'@param data_lats A numeric vector indicating the latitudes of the data. +#'@param data_lons A numeric vector indicating the longitudes of the data. +#'@param mask_sea_land An array with dimensions [lat_dim = data_lats, lon_dim = +#' data_lons] for blending 'data_tas' and 'data_tos'. +#'@param sea_value A numeric value indicating the sea grid points in +#' 'mask_sea_land'. +#'@param type A character string indicating the type of data ('dcpp' for +#' decadal predictions, 'hist' for historical simulations, or 'obs' for +#' observations or reanalyses). +#'@param lat_dim A character string of the name of the latitude dimension. The +#' default value is 'lat'. +#'@param lon_dim A character string of the name of the longitude dimension. The +#' default value is 'lon'. +#'@param mask An array of a mask (with 0's in the grid points that have to be +#' masked) or NULL (i.e., no mask is used). This parameter allows to remove +#' the values over land in case the dataset is a combination of surface air +#' temperature over land and sea surface temperature over the ocean. Also, it +#' can be used to mask those grid points that are missing in the observational +#' dataset for a fair comparison between the forecast system and the reference +#' dataset. The default value is NULL. +#'@param monini An integer indicating the month in which the forecast system is +#' initialized. Only used when parameter 'type' is 'dcpp'. The default value +#' is 11, i.e., initialized in November. +#'@param fmonth_dim A character string indicating the name of the forecast +#' month dimension. Only used if parameter 'type' is 'dcpp'. The default value +#' is 'fmonth'. +#'@param sdate_dim A character string indicating the name of the start date +#' dimension. Only used if parameter 'type' is 'dcpp'. The default value is +#' 'sdate'. +#'@param indices_for_clim A numeric vector of the indices of the years to +#' compute the climatology for calculating the anomalies, or NULL so the +#' climatology is calculated over the whole period. If the data are already +#' 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. +#'@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'. +#'@param month_dim A character string indicating the name of the month +#' dimension. The default value is 'month'. Only used if parameter 'type' is +#' 'hist' or 'obs'. +#'@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'. #' -#'@return The GMST anomalies as function of the sdate, forecast year, and member (in case of decadal predictions); -#' as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +#'@return A numerical array of the GMST anomalies 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). #' #'@examples #' ## Observations or reanalyses -#' obs_tas = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) -#' obs_tos = array(2:101, dim = c(year = 5, lat = 19, lon = 37, month = 12)) -#' mask_sea_land = array(c(1,0,1), dim = c(lat = 19, lon = 37)) -#' sea_value = 1 -#' lat = seq(-90, 90, 10) -#' lon = seq(0, 360, 10) -#' index_obs = GMST(data_tas = obs_tas, data_tos = obs_tos, data_lats = lat, data_lons = lon, type = 'obs', mask_sea_land = mask_sea_land, sea_value = sea_value) +#' obs_tas <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +#' obs_tos <- array(2:101, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +#' mask_sea_land <- array(c(1,0,1), dim = c(lat = 19, lon = 37)) +#' sea_value <- 1 +#' lat <- seq(-90, 90, 10) +#' lon <- seq(0, 360, 10) +#' index_obs <- GMST(data_tas = obs_tas, data_tos = obs_tos, data_lats = lat, data_lons = lon, type = 'obs', mask_sea_land = mask_sea_land, sea_value = sea_value) #' #' ## Historical simulations -#' hist_tas = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) -#' hist_tos = array(2:101, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) -#' mask_sea_land = array(c(1,0,1), dim = c(lat = 19, lon = 37)) -#' sea_value = 1 -#' lat = seq(-90, 90, 10) -#' lon = seq(0, 360, 10) -#' index_hist = GMST(data_tas = hist_tas, data_tos = hist_tos, data_lats = lat, data_lons = lon, type = 'hist', mask_sea_land = mask_sea_land, sea_value = sea_value) +#' hist_tas <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +#' hist_tos <- array(2:101, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +#' mask_sea_land <- array(c(1,0,1), dim = c(lat = 19, lon = 37)) +#' sea_value <- 1 +#' lat <- seq(-90, 90, 10) +#' lon <- seq(0, 360, 10) +#' index_hist <- GMST(data_tas = hist_tas, data_tos = hist_tos, data_lats = lat, data_lons = lon, type = 'hist', mask_sea_land = mask_sea_land, sea_value = sea_value) #' #' ## Decadal predictions -#' dcpp_tas = array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) -#' dcpp_tos = array(2:101, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) -#' mask_sea_land = array(c(1,0,1), dim = c(lat = 19, lon = 37)) -#' sea_value = 1 -#' lat = seq(-90, 90, 10) -#' lon = seq(0, 360, 10) -#' index_dcpp = GMST(data_tas = dcpp_tas, data_tos = dcpp_tos, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1, mask_sea_land = mask_sea_land, sea_value = sea_value) +#' dcpp_tas <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +#' dcpp_tos <- array(2:101, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +#' mask_sea_land <- array(c(1,0,1), dim = c(lat = 19, lon = 37)) +#' sea_value <- 1 +#' lat <- seq(-90, 90, 10) +#' lon <- seq(0, 360, 10) +#' index_dcpp <- GMST(data_tas = dcpp_tas, data_tos = dcpp_tos, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1, mask_sea_land = mask_sea_land, sea_value = sea_value) #' -#'@import ClimProjDiags +#'@author Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} +#'@author Roberto Bilbao, \email{roberto.bilbao@bsc.es} +#'@author Núria Pérez-Zanón, \email{nuria.perez@bsc.es} +#' +#'@importFrom ClimProjDiags WeightedMean #'@import multiApply #'@export GMST <- function(data_tas, data_tos, data_lats, data_lons, mask_sea_land, sea_value, type, mask = NULL, lat_dim = 'lat', lon_dim = 'lon', monini = 11, fmonth_dim = 'fmonth', sdate_dim = 'sdate', indices_for_clim = NULL, - year_dim = 'year', month_dim = 'month', member_dim = 'member'){ + year_dim = 'year', month_dim = 'month', member_dim = 'member') { - ## Checkings - if (!is.array(data_tas)){ - stop('data_tas must be an array') - } - if (!is.array(data_tos)){ - stop('data_tos must be an array') - } - if (!identical(dim(data_tas),dim(data_tos))){ - stop('The dimension of data_tas and data_tos must be identical') + ## Input Checks + # data_tas and data_tos + if (is.null(data_tas) | is.null(data_tos)) { + stop("Parameter 'data_tas' and 'data_tos' cannot be NULL.") } - if(!class(data_lats)=='numeric'){ - stop('data_lats must be a numeric vector') + if (!is.numeric(data_tas) | !is.numeric(data_tos)) { + stop("Parameter 'data_tas' and 'data_tos' must be a numeric array.") } - if (!class(data_lons)=='numeric'){ - stop('data_lons must be a numeric vector or NULL') + if (!identical(dim(data_tas), dim(data_tos))) { + stop("The dimension of data_tas and data_tos must be identical.") } - if (!type %in% c('dcpp','hist','obs')){ - stop("type must be 'dcpp', 'hist' or 'obs'") + # data_lats and data_lons part1 + if (!(class(data_lats) == 'numeric' | class(data_lats) == 'integer')) { + stop("Parameter 'data_lats' must be a numeric vector.") } - if (!is.character(lat_dim)){ - stop('lat_dim must be a string') + if (!(class(data_lons) == 'numeric' | class(data_lons) == 'integer')) { + stop("Parameter 'data_lons' must be a numeric vector.") } - if (!is.character(lon_dim)){ - stop('lon_dim must be a string') + # lat_dim + if (!(is.character(lat_dim) & length(lat_dim) == 1)) { + stop("Parameter 'lat_dim' must be a character string.") } - if (!monini){ - stop("monini must be an integer from 1 to 12") + if (!lat_dim %in% names(dim(data_tas)) | !lat_dim %in% names(dim(data_tos))) { + stop("Parameter 'lat_dim' is not found in 'data_tas' or 'data_tos' dimension.") } - if (!is.character(fmonth_dim)){ - stop('fmonth_dim must be a string') + # lon_dim + if (!(is.character(lon_dim) & length(lon_dim) == 1)) { + stop("Parameter 'lon_dim' must be a character string.") } - if (!is.character(sdate_dim)){ - stop('sdate_dim must be a string') + if (!lon_dim %in% names(dim(data_tas)) | !lon_dim %in% names(dim(data_tos))) { + stop("Parameter 'lon_dim' is not found in 'data_tas' or 'data_tos' dimension.") } - if (!is.null(indices_for_clim) & !class(indices_for_clim) %in% c('numeric','integer') & !is.null(indices_for_clim) & !isFALSE(indices_for_clim)){ - stop("indices_for_clim must be a numeric vector, NULL to compute the anomalies based on the whole period, or FALSE if data are already anomalies") + # data_lats and data_lons part2 + if (dim(data_tas)[lat_dim] != length(data_lats) | + dim(data_tos)[lat_dim] != length(data_lats)) { + stop(paste0("The latitude dimension of parameter 'data_tas' and 'data_tos'", + " must be the same length of parameter 'data_lats'.")) } - if (!is.character(year_dim)){ - stop('year_dim must be a string') - } - if (!is.character(month_dim)){ - stop('month_dim must be a string') + if (dim(data_tas)[lon_dim] != length(data_lons) | + dim(data_tos)[lon_dim] != length(data_lons)) { + stop(paste0("The longitude dimension of parameter 'data_tas' and 'data_tos'", + " must be the same length of parameter 'data_lons'.")) } + # mask_sea_land if (!is.array(mask_sea_land)){ - stop('mask_sea_land must be an array with c(lat_dim,lon_dim) dimensions') - } else if (!identical(names(dim(mask_sea_land)),c(lat_dim,lon_dim))){ - stop('mask_sea_land must be an array with c(lat_dim,lon_dim) dimensions') - } else if (!identical(as.integer(dim(mask_sea_land)),c(length(data_lats),length(data_lons)))){ - stop('mask_sea_land dimension must be equal to the data_lats and data_lons lengths') + stop("Parameter 'mask_sea_land' must be an array with dimensions [lat_dim, lon_dim].") + } else if (!identical(names(dim(mask_sea_land)), c(lat_dim, lon_dim))) { + stop("Parameter 'mask_sea_land' must be an array with dimensions [lat_dim, lon_dim].") + } else if (!identical(as.integer(dim(mask_sea_land)), + c(length(data_lats), length(data_lons)))) { + stop("Parameter 'mask_sea_land' dimensions must be equal to the length of 'data_lats' and 'data_lons'.") } ## combination of tas and tos (data) mask_tas_tos <- function(data_tas, data_tos, mask_sea_land, sea_value) { data <- data_tas - data[mask_sea_land==sea_value] <- data_tos[mask_sea_land==sea_value] + data[mask_sea_land == sea_value] <- data_tos[mask_sea_land == sea_value] return(data) } mask_sea_land <- s2dv::Reorder(data = mask_sea_land, order = c(lat_dim,lon_dim)) - data <- multiApply::Apply(data = list(data_tas, data_tos), target_dims = c(lat_dim,lon_dim), - fun = mask_tas_tos, mask_sea_land = mask_sea_land, sea_value = sea_value)$output1 + data <- multiApply::Apply(data = list(data_tas, data_tos), + target_dims = c(lat_dim, lon_dim), + fun = mask_tas_tos, mask_sea_land = mask_sea_land, + sea_value = sea_value)$output1 data <- drop(data) - rm(data_tas, data_tos) - if (!is.null(mask)){ - if (is.array(mask) & identical(names(dim(mask)),c(lat_dim,lon_dim)) & identical(as.integer(dim(mask)),c(length(data_lats),length(data_lons)))){ + # sea_value + if (!is.numeric(sea_value) | length(sea_value) != 1) { + stop("Parameter 'sea_value' must be a numeric value.") + } + # type + if (!type %in% c('dcpp', 'hist', 'obs')) { + stop("Parameter 'type' must be 'dcpp', 'hist', or 'obs'.") + } + # mask + if (!is.null(mask)) { + if (is.array(mask) & identical(names(dim(mask)), c(lat_dim,lon_dim)) & + identical(as.integer(dim(mask)), c(length(data_lats), length(data_lons)))) { ## To mask those grid point that are missing in the observations - mask <- s2dv::Reorder(data = mask, order = c(lat_dim,lon_dim)) - fun_mask <- function(data, mask){ + mask <- s2dv::Reorder(data = mask, order = c(lat_dim, lon_dim)) + fun_mask <- function(data, mask) { data[mask == 0] <- NA return(data) } - data <- multiApply::Apply(data = data, target_dims = c(lat_dim,lon_dim), fun = fun_mask, mask = mask)$output1 + data <- multiApply::Apply(data = data, target_dims = c(lat_dim, lon_dim), + fun = fun_mask, mask = mask)$output1 } else { - stop('mask must be NULL (no mask) or an array with c(lat_dim,lon_dim) dimensions and 0 in those grid points that have to be masked') + stop(paste0("Parameter 'mask' must be NULL (no mask) or a numerical array ", + "with c(lat_dim, lon_dim) dimensions and 0 in those grid ", + "points that have to be masked.")) } } - - data <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = NULL, - londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) + # monini + if (type == 'dcpp') { + if (!is.numeric(monini) | monini %% 1 != 0 | monini < 1 | + monini > 12) { + stop("Parameter 'monini' must be an integer from 1 to 12.") + } + } + # fmonth_dim + if (type == 'dcpp') { + if (!(is.character(fmonth_dim) & length(fmonth_dim) == 1)) { + stop("Parameter 'fmonth_dim' must be a character string.") + } + if (!fmonth_dim %in% names(dim(data_tas)) | !fmonth_dim %in% names(dim(data_tos))) { + stop("Parameter 'fmonth_dim' is not found in 'data_tas' or 'data_tos' dimension.") + } + } + # sdate_dim + if (type == 'dcpp') { + if (!(is.character(sdate_dim) & length(sdate_dim) == 1)) { + stop("Parameter 'sdate_dim' must be a character string.") + } + if (!sdate_dim %in% names(dim(data_tas)) | !sdate_dim %in% names(dim(data_tos))) { + stop("Parameter 'sdate_dim' is not found in 'data_tas' or 'data_tos' dimension.") + } + } + # indices_for_clim + if (!is.null(indices_for_clim)) { + if (!class(indices_for_clim) %in% c('numeric', 'integer') + & !(is.logical(indices_for_clim) & !any(indices_for_clim))) { + stop(paste0("The parameter 'indices_for_clim' must be a numeric vector ", + "or NULL to compute the anomalies based on the whole period, ", + "or FALSE if data are already anomalies")) + } + } + # year_dim + if (type == 'hist' | type == 'obs') { + if (!(is.character(year_dim) & length(year_dim) == 1)) { + stop("Parameter 'year_dim' must be a character string.") + } + if (!year_dim %in% names(dim(data_tas)) | !year_dim %in% names(dim(data_tos))) { + stop("Parameter 'year_dim' is not found in 'data_tas' or 'data_tos' dimension.") + } + } + # month_dim + if (type == 'hist' | type == 'obs') { + if (!(is.character(month_dim) & length(month_dim) == 1)) { + stop("Parameter 'month_dim' must be a character string.") + } + if (!month_dim %in% names(dim(data_tas)) | !month_dim %in% names(dim(data_tos))) { + stop("Parameter 'month_dim' is not found in 'data_tas' or 'data_tos' 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_tas)) | !member_dim %in% names(dim(data_tos))) { + stop("Parameter 'member_dim' is not found in 'data_tas' or 'data_tos' dimension.") + } + } + + rm(data_tas, data_tos) + + data <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, + region = NULL, + londim = which(names(dim(data)) == lon_dim), + latdim = which(names(dim(data)) == lat_dim)) INDEX <- .Indices(data = data, type = type, monini = monini, indices_for_clim = indices_for_clim, fmonth_dim = fmonth_dim, diff --git a/R/GSAT.R b/R/GSAT.R index ced384c..b9d0f6f 100644 --- a/R/GSAT.R +++ b/R/GSAT.R @@ -1,113 +1,218 @@ #'Compute the Global Surface Air Temperature (GSAT) anomalies #' -#'@description The Global Surface Air Temperature (GSAT) anomalies are computed as the +#'The Global Surface Air Temperature (GSAT) anomalies are computed as the #'weighted-averaged surface air temperature anomalies over the global region. -#'@author Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} -#'@author Núria Pérez-Zanón, \email{nuria.perez@bsc.es} -#'@author Roberto Bilbao, \email{roberto.bilbao@bsc.es} #' -#'@param data Data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), -#' with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). -#' This data has to be provided, at least, over the whole region needed to compute the index. -#'@param data_lats An array with the latitudes of the data. -#'@param data_lons An array with the longitudes of the data. -#'@param type A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses). -#'@param lat_dim A string with the name of the latitude dimension ('lat' by default). -#'@param lon_dim A string with the name of the longitude dimension ('lon' by default). -#'@param mask An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a -#' combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the -#' observational dataset for a fair comparison between the forecast system and the reference dataset. -#'@param monini Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'. -#'@param fmonth_dim A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'. -#'@param sdate_dim A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'. -#'@param indices_for_clim Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). -#' In case of type='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. -#'@param year_dim A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'. -#'@param month_dim A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'. -#'@param member_dim A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'. +#'@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 +#' 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. +#'@param type A character string indicating the type of data ('dcpp' for +#' decadal predictions, 'hist' for historical simulations, or 'obs' for +#' observations or reanalyses). +#'@param lat_dim A character string of the name of the latitude dimension. The +#' default value is 'lat'. +#'@param lon_dim A character string of the name of the longitude dimension. The +#' default value is 'lon'. +#'@param mask An array of a mask (with 0's in the grid points that have to be +#' masked) or NULL (i.e., no mask is used). This parameter allows to remove +#' the values over land in case the dataset is a combination of surface air +#' temperature over land and sea surface temperature over the ocean. Also, it +#' can be used to mask those grid points that are missing in the observational +#' dataset for a fair comparison between the forecast system and the reference +#' dataset. The default value is NULL. +#'@param monini An integer indicating the month in which the forecast system is +#' initialized. Only used when parameter 'type' is 'dcpp'. The default value +#' is 11, i.e., initialized in November. +#'@param fmonth_dim A character string indicating the name of the forecast +#' month dimension. Only used if parameter 'type' is 'dcpp'. The default value +#' is 'fmonth'. +#'@param sdate_dim A character string indicating the name of the start date +#' dimension. Only used if parameter 'type' is 'dcpp'. The default value is +#' 'sdate'. +#'@param indices_for_clim A numeric vector of the indices of the years to +#' compute the climatology for calculating the anomalies, or NULL so the +#' climatology is calculated over the whole period. If the data are already +#' 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. +#'@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'. +#'@param month_dim A character string indicating the name of the month +#' dimension. The default value is 'month'. Only used if parameter 'type' is +#' 'hist' or 'obs'. +#'@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'. #' -#'@return The GSAT anomalies as function of the sdate, forecast year, and member (in case of decadal predictions); -#' as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +#'@return A numerical array of the GSAT anomalies 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). #' #'@examples #' ## Observations or reanalyses -#' obs = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) -#' lat = seq(-90, 90, 10) -#' lon = seq(0, 360, 10) -#' index_obs = GSAT(data = obs, data_lats = lat, data_lons = lon, type = 'obs') +#' obs <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +#' lat <- seq(-90, 90, 10) +#' lon <- seq(0, 360, 10) +#' index_obs <- GSAT(data = obs, data_lats = lat, data_lons = lon, type = 'obs') #' #' ## Historical simulations -#' hist = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) -#' lat = seq(-90, 90, 10) -#' lon = seq(0, 360, 10) -#' index_hist = GSAT(data = hist, data_lats = lat, data_lons = lon, type = 'hist') +#' hist <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +#' lat <- seq(-90, 90, 10) +#' lon <- seq(0, 360, 10) +#' index_hist <- GSAT(data = hist, data_lats = lat, data_lons = lon, type = 'hist') #' #' ## Decadal predictions -#' dcpp = array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) -#' lat = seq(-90, 90, 10) -#' lon = seq(0, 360, 10) -#' index_dcpp = GSAT(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) +#' dcpp <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +#' lat <- seq(-90, 90, 10) +#' lon <- seq(0, 360, 10) +#' index_dcpp <- GSAT(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) +#' +#'@author Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} +#'@author Núria Pérez-Zanón, \email{nuria.perez@bsc.es} +#'@author Roberto Bilbao, \email{roberto.bilbao@bsc.es} #' -#'@import ClimProjDiags +#'@importFrom ClimProjDiags WeightedMean #'@import multiApply #'@export GSAT <- 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'){ + indices_for_clim = NULL, year_dim = 'year', month_dim = 'month', + member_dim = 'member') { - ## Checkings - if (!is.array(data)){ - stop('data must be an array') + ## Input Checks + # data + if (is.null(data)) { + stop("Parameter 'data' cannot be NULL.") } - if(!class(data_lats)=='numeric'){ - stop('data_lats must be a numeric vector') + if (!is.numeric(data)) { + stop("Parameter 'data' must be a numeric array.") } - if (!class(data_lons)=='numeric'){ - stop('data_lons must be a numeric vector or NULL') + # data_lats and data_lons part1 + if (!(class(data_lats) == 'numeric' | class(data_lats) == 'integer')) { + stop("Parameter 'data_lats' must be a numeric vector.") } - if (!type %in% c('dcpp','hist','obs')){ - stop("type must be 'dcpp', 'hist' or 'obs'") + if (!(class(data_lons) == 'numeric' | class(data_lons) == 'integer')) { + stop("Parameter 'data_lons' must be a numeric vector.") } - if (!is.character(lat_dim)){ - stop('lat_dim must be a string') + # type + if (!type %in% c('dcpp', 'hist', 'obs')) { + stop("Parameter 'type' must be 'dcpp', 'hist', or 'obs'.") } - if (!is.character(lon_dim)){ - stop('lon_dim must be a string') + # lat_dim + if (!(is.character(lat_dim) & length(lat_dim) == 1)) { + stop("Parameter 'lat_dim' must be a character string.") } - if (!monini){ - stop("monini must be an integer from 1 to 12") + if (!lat_dim %in% names(dim(data))) { + stop("Parameter 'lat_dim' is not found in 'data' dimension.") } - if (!is.character(fmonth_dim)){ - stop('fmonth_dim must be a string') + # lon_dim + if (!(is.character(lon_dim) & length(lon_dim) == 1)) { + stop("Parameter 'lon_dim' must be a character string.") } - if (!is.character(sdate_dim)){ - stop('sdate_dim must be a string') + if (!lon_dim %in% names(dim(data))) { + stop("Parameter 'lon_dim' is not found in 'data' dimension.") } - if (!is.null(indices_for_clim) & !class(indices_for_clim) %in% c('numeric','integer') & !is.null(indices_for_clim) & !isFALSE(indices_for_clim)){ - stop("indices_for_clim must be a numeric vector, NULL to compute the anomalies based on the whole period, or FALSE if data are already anomalies") + # data_lats and data_lons part2 + if (dim(data)[lat_dim] != length(data_lats)){ + stop(paste0("The latitude dimension of parameter 'data' must be the same", + " length of parameter 'data_lats'.")) } - if (!is.character(year_dim)){ - stop('year_dim must be a string') + if (dim(data)[lon_dim] != length(data_lons)){ + stop(paste0("The longitude dimension of parameter 'data' must be the same", + " length of parameter 'data_lons'.")) } - if (!is.character(month_dim)){ - stop('month_dim must be a string') - } - if (!is.null(mask)){ - if (is.array(mask) & identical(names(dim(mask)),c(lat_dim,lon_dim)) & identical(as.integer(dim(mask)),c(length(data_lats),length(data_lons)))){ + # mask + if (!is.null(mask)) { + if (is.array(mask) & identical(names(dim(mask)), c(lat_dim,lon_dim)) & + identical(as.integer(dim(mask)), c(length(data_lats), length(data_lons)))) { ## To mask those grid point that are missing in the observations - mask <- s2dv::Reorder(data = mask, order = c(lat_dim,lon_dim)) + mask <- s2dv::Reorder(data = mask, order = c(lat_dim, lon_dim)) fun_mask <- function(data, mask){ data[mask == 0] <- NA return(data) } - data <- multiApply::Apply(data = data, target_dims = c(lat_dim,lon_dim), fun = fun_mask, mask = mask)$output1 + data <- multiApply::Apply(data = data, target_dims = c(lat_dim, lon_dim), + fun = fun_mask, mask = mask)$output1 } else { - stop('mask must be NULL (no mask) or an array with c(lat_dim,lon_dim) dimensions and 0 in those grid points that have to be masked') + stop(paste0("Parameter 'mask' must be NULL (no mask) or a numerical array ", + "with c(lat_dim, lon_dim) dimensions and 0 in those grid ", + "points that have to be masked.")) } } - - data <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = NULL, - londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) + # monini + if (type == 'dcpp') { + if (!is.numeric(monini) | monini %% 1 != 0 | monini < 1 | + monini > 12) { + stop("Parameter 'monini' must be an integer from 1 to 12.") + } + } + # fmonth_dim + if (type == 'dcpp') { + if (!(is.character(fmonth_dim) & length(fmonth_dim) == 1)) { + stop("Parameter 'fmonth_dim' must be a character string.") + } + if (!fmonth_dim %in% names(dim(data))) { + stop("Parameter 'fmonth_dim' is not found in 'data' dimension.") + } + } + # sdate_dim + if (type == 'dcpp') { + if (!(is.character(sdate_dim) & length(sdate_dim) == 1)) { + stop("Parameter 'sdate_dim' must be a character string.") + } + if (!sdate_dim %in% names(dim(data))) { + stop("Parameter 'sdate_dim' is not found in 'data' dimension.") + } + } + # indices_for_clim + if (!is.null(indices_for_clim)) { + if (!class(indices_for_clim) %in% c('numeric', 'integer') + & !(is.logical(indices_for_clim) & !any(indices_for_clim))) { + stop(paste0("The parameter 'indices_for_clim' must be a numeric vector ", + "or NULL to compute the anomalies based on the whole period, ", + "or FALSE if data are already anomalies")) + } + } + # year_dim + if (type == 'hist' | type == 'obs') { + if (!(is.character(year_dim) & length(year_dim) == 1)) { + stop("Parameter 'year_dim' must be a character string.") + } + if (!year_dim %in% names(dim(data))) { + stop("Parameter 'year_dim' is not found in 'data' dimension.") + } + } + # month_dim + if (type == 'hist' | type == 'obs') { + if (!(is.character(month_dim) & length(month_dim) == 1)) { + stop("Parameter 'month_dim' must be a character string.") + } + if (!month_dim %in% names(dim(data))) { + 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.") + } + } + + data <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, + region = NULL, + londim = which(names(dim(data)) == lon_dim), + latdim = which(names(dim(data)) == lat_dim)) INDEX <- .Indices(data = data, type = type, monini = monini, indices_for_clim = indices_for_clim, fmonth_dim = fmonth_dim, diff --git a/R/SPOD.R b/R/SPOD.R index 688a527..d939df7 100644 --- a/R/SPOD.R +++ b/R/SPOD.R @@ -1,127 +1,237 @@ #'Compute the South Pacific Ocean Dipole (SPOD) index #' -#'@description The South Pacific Ocean Dipole (SPOD) index is related to the El Nino-Southern Oscillation (ENSO) and the Inderdecadal Pacific Oscillation (IPO). -#'The SPOD index is computed as the difference of weighted-averaged SST anomalies over 20ºS-48ºS, 165ºE-190ºE (NW pole) and the weighted-averaged SST anomalies over -#'44ºS-65ºS, 220ºE-260ºE (SE pole) (Saurral et al., 2020). -#'@author Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} -#'@author Núria Pérez-Zanón, \email{nuria.perez@bsc.es} -#'@author Roberto Bilbao, \email{roberto.bilbao@bsc.es} +#'The South Pacific Ocean Dipole (SPOD) index is related to the El +#'Nino-Southern Oscillation (ENSO) and the Inderdecadal Pacific Oscillation +#'(IPO). The SPOD index is computed as the difference of weighted-averaged SST +#'anomalies over 20ºS-48ºS, 165ºE-190ºE (NW pole) and the weighted-averaged SST +#' anomalies over 44ºS-65ºS, 220ºE-260ºE (SE pole) (Saurral et al., 2020). #' -#'@param data Data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), -#' with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). -#' This data has to be provided, at least, over the whole region needed to compute the index. -#'@param data_lats An array with the latitudes of the data. -#'@param data_lons An array with the longitudes of the data. -#'@param type A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses). -#'@param lat_dim A string with the name of the latitude dimension ('lat' by default). -#'@param lon_dim A string with the name of the longitude dimension ('lon' by default). -#'@param mask An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a -#' combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the -#' observational dataset for a fair comparison between the forecast system and the reference dataset. -#'@param monini Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'. -#'@param fmonth_dim A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'. -#'@param sdate_dim A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'. -#'@param indices_for_clim Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). -#' In case of type='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. -#'@param year_dim A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'. -#'@param month_dim A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'. -#'@param member_dim A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'. +#'@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 +#' 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. +#'@param type A character string indicating the type of data ('dcpp' for +#' decadal predictions, 'hist' for historical simulations, or 'obs' for +#' observations or reanalyses). +#'@param lat_dim A character string of the name of the latitude dimension. The +#' default value is 'lat'. +#'@param lon_dim A character string of the name of the longitude dimension. The +#' default value is 'lon'. +#'@param mask An array of a mask (with 0's in the grid points that have to be +#' masked) or NULL (i.e., no mask is used). This parameter allows to remove +#' the values over land in case the dataset is a combination of surface air +#' temperature over land and sea surface temperature over the ocean. Also, it +#' can be used to mask those grid points that are missing in the observational +#' dataset for a fair comparison between the forecast system and the reference +#' dataset. The default value is NULL. +#'@param monini An integer indicating the month in which the forecast system is +#' initialized. Only used when parameter 'type' is 'dcpp'. The default value +#' is 11, i.e., initialized in November. +#'@param fmonth_dim A character string indicating the name of the forecast +#' month dimension. Only used if parameter 'type' is 'dcpp'. The default value +#' is 'fmonth'. +#'@param sdate_dim A character string indicating the name of the start date +#' dimension. Only used if parameter 'type' is 'dcpp'. The default value is +#' 'sdate'. +#'@param indices_for_clim A numeric vector of the indices of the years to +#' compute the climatology for calculating the anomalies, or NULL so the +#' climatology is calculated over the whole period. If the data are already +#' 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. +#'@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'. +#'@param month_dim A character string indicating the name of the month +#' dimension. The default value is 'month'. Only used if parameter 'type' is +#' 'hist' or 'obs'. +#'@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'. #' -#'@return The SPOD index as function of the sdate, forecast year, and member (in case of decadal predictions); -#' as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +#'@return A numerical array of the SPOD 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). #' #'@examples #' ## Observations or reanalyses -#' obs = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) -#' lat = seq(-90, 90, 10) -#' lon = seq(0, 360, 10) -#' index_obs = SPOD(data = obs, data_lats = lat, data_lons = lon, type = 'obs') +#' obs <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +#' lat <- seq(-90, 90, 10) +#' lon <- seq(0, 360, 10) +#' index_obs <- SPOD(data = obs, data_lats = lat, data_lons = lon, type = 'obs') #' #' ## Historical simulations -#' hist = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) -#' lat = seq(-90, 90, 10) -#' lon = seq(0, 360, 10) -#' index_hist = SPOD(data = hist, data_lats = lat, data_lons = lon, type = 'hist') +#' hist <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +#' lat <- seq(-90, 90, 10) +#' lon <- seq(0, 360, 10) +#' index_hist <- SPOD(data = hist, data_lats = lat, data_lons = lon, type = 'hist') #' #' ## Decadal predictions -#' dcpp = array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) -#' lat = seq(-90, 90, 10) -#' lon = seq(0, 360, 10) -#' index_dcpp = SPOD(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) +#' dcpp <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +#' lat <- seq(-90, 90, 10) +#' lon <- seq(0, 360, 10) +#' index_dcpp <- SPOD(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) +#' +#'@author Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} +#'@author Núria Pérez-Zanón, \email{nuria.perez@bsc.es} +#'@author Roberto Bilbao, \email{roberto.bilbao@bsc.es} #' -#'@import ClimProjDiags +#'@importFrom ClimProjDiags WeightedMean CombineIndices #'@import multiApply #'@export SPOD <- 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'){ + indices_for_clim = NULL, year_dim = 'year', month_dim = 'month', + member_dim = 'member') { - ## Checkings - if (!is.array(data)){ - stop('data must be an array') + ## Input Checks + # data + if (is.null(data)) { + stop("Parameter 'data' cannot be NULL.") } - if(!class(data_lats)=='numeric'){ - stop('data_lats must be a numeric vector') + if (!is.numeric(data)) { + stop("Parameter 'data' must be a numeric array.") } - if (!class(data_lons)=='numeric'){ - stop('data_lons must be a numeric vector or NULL') + # data_lats and data_lons part1 + if (!(class(data_lats) == 'numeric' | class(data_lats) == 'integer')) { + stop("Parameter 'data_lats' must be a numeric vector.") } - if (!type %in% c('dcpp','hist','obs')){ - stop("type must be 'dcpp', 'hist' or 'obs'") + if (!(class(data_lons) == 'numeric' | class(data_lons) == 'integer')) { + stop("Parameter 'data_lons' must be a numeric vector.") } - if (!is.character(lat_dim)){ - stop('lat_dim must be a string') + # type + if (!type %in% c('dcpp', 'hist', 'obs')) { + stop("Parameter 'type' must be 'dcpp', 'hist', or 'obs'.") } - if (!is.character(lon_dim)){ - stop('lon_dim must be a string') + # lat_dim + if (!(is.character(lat_dim) & length(lat_dim) == 1)) { + stop("Parameter 'lat_dim' must be a character string.") } - if (!monini){ - stop("monini must be an integer from 1 to 12") + if (!lat_dim %in% names(dim(data))) { + stop("Parameter 'lat_dim' is not found in 'data' dimension.") } - if (!is.character(fmonth_dim)){ - stop('fmonth_dim must be a string') + # lon_dim + if (!(is.character(lon_dim) & length(lon_dim) == 1)) { + stop("Parameter 'lon_dim' must be a character string.") } - if (!is.character(sdate_dim)){ - stop('sdate_dim must be a string') + if (!lon_dim %in% names(dim(data))) { + stop("Parameter 'lon_dim' is not found in 'data' dimension.") } - if (!is.null(indices_for_clim) & !class(indices_for_clim) %in% c('numeric','integer') & !is.null(indices_for_clim) & !isFALSE(indices_for_clim)){ - stop("indices_for_clim must be a numeric vector, NULL to compute the anomalies based on the whole period, or FALSE if data are already anomalies") + # data_lats and data_lons part2 + if (dim(data)[lat_dim] != length(data_lats)){ + stop(paste0("The latitude dimension of parameter 'data' must be the same", + " length of parameter 'data_lats'.")) } - if (!is.character(year_dim)){ - stop('year_dim must be a string') + if (dim(data)[lon_dim] != length(data_lons)){ + stop(paste0("The longitude dimension of parameter 'data' must be the same", + " length of parameter 'data_lons'.")) } - if (!is.character(month_dim)){ - stop('month_dim must be a string') - } - if (!is.null(mask)){ - if (is.array(mask) & identical(names(dim(mask)),c(lat_dim,lon_dim)) & identical(as.integer(dim(mask)),c(length(data_lats),length(data_lons)))){ + # mask + if (!is.null(mask)) { + if (is.array(mask) & identical(names(dim(mask)), c(lat_dim,lon_dim)) & + identical(as.integer(dim(mask)), c(length(data_lats), length(data_lons)))) { ## To mask those grid point that are missing in the observations - mask <- s2dv::Reorder(data = mask, order = c(lat_dim,lon_dim)) + mask <- s2dv::Reorder(data = mask, order = c(lat_dim, lon_dim)) fun_mask <- function(data, mask){ data[mask == 0] <- NA return(data) } - data <- multiApply::Apply(data = data, target_dims = c(lat_dim,lon_dim), fun = fun_mask, mask = mask)$output1 + data <- multiApply::Apply(data = data, target_dims = c(lat_dim, lon_dim), + fun = fun_mask, mask = mask)$output1 } else { - stop('mask must be NULL (no mask) or an array with c(lat_dim,lon_dim) dimensions and 0 in those grid points that have to be masked') + stop(paste0("Parameter 'mask' must be NULL (no mask) or a numerical array ", + "with c(lat_dim, lon_dim) dimensions and 0 in those grid ", + "points that have to be masked.")) } } - + # monini + if (type == 'dcpp') { + if (!is.numeric(monini) | monini %% 1 != 0 | monini < 1 | + monini > 12) { + stop("Parameter 'monini' must be an integer from 1 to 12.") + } + } + # fmonth_dim + if (type == 'dcpp') { + if (!(is.character(fmonth_dim) & length(fmonth_dim) == 1)) { + stop("Parameter 'fmonth_dim' must be a character string.") + } + if (!fmonth_dim %in% names(dim(data))) { + stop("Parameter 'fmonth_dim' is not found in 'data' dimension.") + } + } + # sdate_dim + if (type == 'dcpp') { + if (!(is.character(sdate_dim) & length(sdate_dim) == 1)) { + stop("Parameter 'sdate_dim' must be a character string.") + } + if (!sdate_dim %in% names(dim(data))) { + stop("Parameter 'sdate_dim' is not found in 'data' dimension.") + } + } + # indices_for_clim + if (!is.null(indices_for_clim)) { + if (!class(indices_for_clim) %in% c('numeric', 'integer') + & !(is.logical(indices_for_clim) & !any(indices_for_clim))) { + stop(paste0("The parameter 'indices_for_clim' must be a numeric vector ", + "or NULL to compute the anomalies based on the whole period, ", + "or FALSE if data are already anomalies")) + } + } + # year_dim + if (type == 'hist' | type == 'obs') { + if (!(is.character(year_dim) & length(year_dim) == 1)) { + stop("Parameter 'year_dim' must be a character string.") + } + if (!year_dim %in% names(dim(data))) { + stop("Parameter 'year_dim' is not found in 'data' dimension.") + } + } + # month_dim + if (type == 'hist' | type == 'obs') { + if (!(is.character(month_dim) & length(month_dim) == 1)) { + stop("Parameter 'month_dim' must be a character string.") + } + if (!month_dim %in% names(dim(data))) { + 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.") + } + } + ## Regions for IPO_SPOD (Saurral et al., 2020) - lat_min_1 <- -48; lat_max_1 = -20 - lon_min_1 <- 165; lon_max_1 = 190 - lat_min_2 <- -65; lat_max_2 = -44 - lon_min_2 <- 220; lon_max_2 = 260 - regions = NULL + lat_min_1 <- -48; lat_max_1 <- -20 + lon_min_1 <- 165; lon_max_1 <- 190 + lat_min_2 <- -65; lat_max_2 <- -44 + lon_min_2 <- 220; lon_max_2 <- 260 + regions <- NULL regions$reg1 <- c(lon_min_1, lon_max_1, lat_min_1, lat_max_1) regions$reg2 <- c(lon_min_2, lon_max_2, lat_min_2, lat_max_2) - mean_1 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = regions$reg1, - londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) - mean_2 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = regions$reg2, - londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) + mean_1 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, + region = regions$reg1, + londim = which(names(dim(data)) == lon_dim), + latdim = which(names(dim(data)) == lat_dim)) + mean_2 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, + region = regions$reg2, + londim = which(names(dim(data)) == lon_dim), + latdim = which(names(dim(data)) == lat_dim)) - data <- ClimProjDiags::CombineIndices(indices = list(mean_1,mean_2), weights = NULL, operation = 'subtract') # (mean_1 - mean_2) + 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, diff --git a/R/TPI.R b/R/TPI.R index fbf6e59..6733bdd 100644 --- a/R/TPI.R +++ b/R/TPI.R @@ -1,33 +1,62 @@ #'Compute the Tripole Index (TPI) for the Interdecadal Pacific Oscillation (IPO) #' -#'@description The Tripole Index (TPI) for the Interdecadal Pacific Oscillation (IPO) is computed as the difference of weighted-averaged SST anomalies over -#'10ºS-10ºN, 170ºE-270ºE minus the mean of the weighted-averaged SST anomalies over 25ºN-45ºN, 140ºE-215ºE and 50ºS-15ºS, 150ºE-200ºE (Henley et al., 2015). -#'@author Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} -#'@author Núria Pérez-Zanón, \email{nuria.perez@bsc.es} -#'@author Roberto Bilbao, \email{roberto.bilbao@bsc.es} +#'The Tripole Index (TPI) for the Interdecadal Pacific Oscillation (IPO) is +#'computed as the difference of weighted-averaged SST anomalies over 10ºS-10ºN, +#'170ºE-270ºE minus the mean of the weighted-averaged SST anomalies over +#'25ºN-45ºN, 140ºE-215ºE and 50ºS-15ºS, 150ºE-200ºE (Henley et al., 2015). #' -#'@param data Data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), -#' with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). -#' This data has to be provided, at least, over the whole region needed to compute the index. -#'@param data_lats An array with the latitudes of the data. -#'@param data_lons An array with the longitudes of the data. -#'@param type A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses). -#'@param lat_dim A string with the name of the latitude dimension ('lat' by default). -#'@param lon_dim A string with the name of the longitude dimension ('lon' by default). -#'@param mask An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a -#' combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the -#' observational dataset for a fair comparison between the forecast system and the reference dataset. -#'@param monini Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'. -#'@param fmonth_dim A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'. -#'@param sdate_dim A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'. -#'@param indices_for_clim Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). -#' In case of type='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. -#'@param year_dim A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'. -#'@param month_dim A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'. -#'@param member_dim A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'. +#'@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 +#' 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. +#'@param A character string indicating the type of data ('dcpp' for +#' decadal predictions, 'hist' for historical simulations, or 'obs' for +#' observations or reanalyses). +#'@param lat_dim A character string of the name of the latitude dimension. The +#' default value is 'lat'. +#'@param lon_dim A character string of the name of the longitude dimension. The +#' default value is 'lon'. +#'@param mask An array of a mask (with 0's in the grid points that have to be +#' masked) or NULL (i.e., no mask is used). This parameter allows to remove +#' the values over land in case the dataset is a combination of surface air +#' temperature over land and sea surface temperature over the ocean. Also, it +#' can be used to mask those grid points that are missing in the observational +#' dataset for a fair comparison between the forecast system and the reference +#' dataset. The default value is NULL. +#'@param monini An integer indicating the month in which the forecast system is +#' initialized. Only used when parameter 'type' is 'dcpp'. The default value +#' is 11, i.e., initialized in November. +#'@param fmonth_dim A character string indicating the name of the forecast +#' month dimension. Only used if parameter 'type' is 'dcpp'. The default value +#' is 'fmonth'. +#'@param sdate_dim A character string indicating the name of the start date +#' dimension. Only used if parameter 'type' is 'dcpp'. The default value is +#' 'sdate'. +#'@param indices_for_clim A numeric vector of the indices of the years to +#' compute the climatology for calculating the anomalies, or NULL so the +#' climatology is calculated over the whole period. If the data are already +#' 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. +#'@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'. +#'@param month_dim A character string indicating the name of the month +#' dimension. The default value is 'month'. Only used if parameter 'type' is +#' 'hist' or 'obs'. +#'@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'. #' -#'@return The TPI as function of the sdate, forecast year, and member (in case of decadal predictions); -#' as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +#'@return A numerical array of the TPI 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). #' #'@examples #' ## Observations or reanalyses @@ -47,62 +76,138 @@ #' lat = seq(-90, 90, 10) #' lon = seq(0, 360, 10) #' index_dcpp = TPI(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) +#' +#'@author Carlos Delgado-Torres, \email{carlos.delgado@bsc.es} +#'@author Núria Pérez-Zanón, \email{nuria.perez@bsc.es} +#'@author Roberto Bilbao, \email{roberto.bilbao@bsc.es} #' -#'@import ClimProjDiags +#'@importFrom ClimProjDiags WeightedMean CombineIndices #'@import multiApply #'@export TPI <- function(data, data_lats, data_lons, index, 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'){ + indices_for_clim = NULL, year_dim = 'year', month_dim = 'month', + member_dim = 'member') { - ## Checkings - if (!is.array(data)){ - stop('data must be an array') + ## Input Checks + # data + if (is.null(data)) { + stop("Parameter 'data' cannot be NULL.") } - if(!class(data_lats)=='numeric'){ - stop('data_lats must be a numeric vector') + if (!is.numeric(data)) { + stop("Parameter 'data' must be a numeric array.") } - if (!class(data_lons)=='numeric'){ - stop('data_lons must be a numeric vector or NULL') + # data_lats and data_lons part1 + if (!(class(data_lats) == 'numeric' | class(data_lats) == 'integer')) { + stop("Parameter 'data_lats' must be a numeric vector.") } - if (!type %in% c('dcpp','hist','obs')){ - stop("type must be 'dcpp', 'hist' or 'obs'") + if (!(class(data_lons) == 'numeric' | class(data_lons) == 'integer')) { + stop("Parameter 'data_lons' must be a numeric vector.") } - if (!is.character(lat_dim)){ - stop('lat_dim must be a string') + # type + if (!type %in% c('dcpp', 'hist', 'obs')) { + stop("Parameter 'type' must be 'dcpp', 'hist', or 'obs'.") } - if (!is.character(lon_dim)){ - stop('lon_dim must be a string') + # lat_dim + if (!(is.character(lat_dim) & length(lat_dim) == 1)) { + stop("Parameter 'lat_dim' must be a character string.") } - if (!monini){ - stop("monini must be an integer from 1 to 12") + if (!lat_dim %in% names(dim(data))) { + stop("Parameter 'lat_dim' is not found in 'data' dimension.") } - if (!is.character(fmonth_dim)){ - stop('fmonth_dim must be a string') + # lon_dim + if (!(is.character(lon_dim) & length(lon_dim) == 1)) { + stop("Parameter 'lon_dim' must be a character string.") } - if (!is.character(sdate_dim)){ - stop('sdate_dim must be a string') + if (!lon_dim %in% names(dim(data))) { + stop("Parameter 'lon_dim' is not found in 'data' dimension.") } - if (!is.null(indices_for_clim) & !class(indices_for_clim) %in% c('numeric','integer') & !is.null(indices_for_clim) & !isFALSE(indices_for_clim)){ - stop("indices_for_clim must be a numeric vector, NULL to compute the anomalies based on the whole period, or FALSE if data are already anomalies") + # data_lats and data_lons part2 + if (dim(data)[lat_dim] != length(data_lats)){ + stop(paste0("The latitude dimension of parameter 'data' must be the same", + " length of parameter 'data_lats'.")) } - if (!is.character(year_dim)){ - stop('year_dim must be a string') + if (dim(data)[lon_dim] != length(data_lons)){ + stop(paste0("The longitude dimension of parameter 'data' must be the same", + " length of parameter 'data_lons'.")) } - if (!is.character(month_dim)){ - stop('month_dim must be a string') - } - if (!is.null(mask)){ - if (is.array(mask) & identical(names(dim(mask)),c(lat_dim,lon_dim)) & identical(as.integer(dim(mask)),c(length(data_lats),length(data_lons)))){ + # mask + if (!is.null(mask)) { + if (is.array(mask) & identical(names(dim(mask)), c(lat_dim,lon_dim)) & + identical(as.integer(dim(mask)), c(length(data_lats), length(data_lons)))) { ## To mask those grid point that are missing in the observations - mask <- s2dv::Reorder(data = mask, order = c(lat_dim,lon_dim)) + mask <- s2dv::Reorder(data = mask, order = c(lat_dim, lon_dim)) fun_mask <- function(data, mask){ data[mask == 0] <- NA return(data) } - data <- multiApply::Apply(data = data, target_dims = c(lat_dim,lon_dim), fun = fun_mask, mask = mask)$output1 + data <- multiApply::Apply(data = data, target_dims = c(lat_dim, lon_dim), + fun = fun_mask, mask = mask)$output1 } else { - stop('mask must be NULL (no mask) or an array with c(lat_dim,lon_dim) dimensions and 0 in those grid points that have to be masked') + stop(paste0("Parameter 'mask' must be NULL (no mask) or a numerical array ", + "with c(lat_dim, lon_dim) dimensions and 0 in those grid ", + "points that have to be masked.")) + } + } + # monini + if (type == 'dcpp') { + if (!is.numeric(monini) | monini %% 1 != 0 | monini < 1 | + monini > 12) { + stop("Parameter 'monini' must be an integer from 1 to 12.") + } + } + # fmonth_dim + if (type == 'dcpp') { + if (!(is.character(fmonth_dim) & length(fmonth_dim) == 1)) { + stop("Parameter 'fmonth_dim' must be a character string.") + } + if (!fmonth_dim %in% names(dim(data))) { + stop("Parameter 'fmonth_dim' is not found in 'data' dimension.") + } + } + # sdate_dim + if (type == 'dcpp') { + if (!(is.character(sdate_dim) & length(sdate_dim) == 1)) { + stop("Parameter 'sdate_dim' must be a character string.") + } + if (!sdate_dim %in% names(dim(data))) { + stop("Parameter 'sdate_dim' is not found in 'data' dimension.") + } + } + # indices_for_clim + if (!is.null(indices_for_clim)) { + if (!class(indices_for_clim) %in% c('numeric', 'integer') + & !(is.logical(indices_for_clim) & !any(indices_for_clim))) { + stop(paste0("The parameter 'indices_for_clim' must be a numeric vector ", + "or NULL to compute the anomalies based on the whole period, ", + "or FALSE if data are already anomalies")) + } + } + # year_dim + if (type == 'hist' | type == 'obs') { + if (!(is.character(year_dim) & length(year_dim) == 1)) { + stop("Parameter 'year_dim' must be a character string.") + } + if (!year_dim %in% names(dim(data))) { + stop("Parameter 'year_dim' is not found in 'data' dimension.") + } + } + # month_dim + if (type == 'hist' | type == 'obs') { + if (!(is.character(month_dim) & length(month_dim) == 1)) { + stop("Parameter 'month_dim' must be a character string.") + } + if (!month_dim %in% names(dim(data))) { + 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.") } } @@ -113,21 +218,29 @@ TPI <- function(data, data_lats, data_lons, index, type, lat_dim = 'lat', lon_di lon_min_2 <- 170; lon_max_2 <- 270 lat_min_3 <- -50; lat_max_3 <- -15 lon_min_3 <- 150; lon_max_3 <- 200 - regions = NULL + regions <- NULL regions$reg1 <- c(lon_min_1, lon_max_1, lat_min_1, lat_max_1) regions$reg2 <- c(lon_min_2, lon_max_2, lat_min_2, lat_max_2) regions$reg3 <- c(lon_min_3, lon_max_3, lat_min_3, lat_max_3) - mean_1 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = regions$reg1, - londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) - mean_2 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = regions$reg2, - londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) - mean_3 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = regions$reg3, - londim = which(names(dim(data))==lon_dim), latdim = which(names(dim(data))==lat_dim)) + mean_1 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, + region = regions$reg1, + londim = which(names(dim(data)) == lon_dim), + latdim = which(names(dim(data)) == lat_dim)) + mean_2 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, + region = regions$reg2, + londim = which(names(dim(data)) == lon_dim), + latdim = which(names(dim(data)) == lat_dim)) + mean_3 <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, + region = regions$reg3, + londim = which(names(dim(data)) == lon_dim), + latdim = which(names(dim(data)) == lat_dim)) - mean_1_3 <- ClimProjDiags::CombineIndices(indices = list(mean_1,mean_3), weights = NULL, operation = 'mean') + mean_1_3 <- ClimProjDiags::CombineIndices(indices = list(mean_1, mean_3), + weights = NULL, operation = 'mean') - data <- ClimProjDiags::CombineIndices(indices = list(mean_2, mean_1_3), weights = NULL, operation = 'subtract') # mean_2 - ((mean_1 + mean_3)/2) + data <- ClimProjDiags::CombineIndices(indices = list(mean_2, mean_1_3), + weights = NULL, operation = 'subtract') # mean_2 - ((mean_1 + mean_3)/2) INDEX <- .Indices(data = data, type = type, monini = monini, indices_for_clim = indices_for_clim, fmonth_dim = fmonth_dim, diff --git a/man/AMV.Rd b/man/AMV.Rd index 3e87dc7..7cb7503 100644 --- a/man/AMV.Rd +++ b/man/AMV.Rd @@ -10,13 +10,12 @@ AMV(data, data_lats, data_lons, type, lat_dim = "lat", lon_dim = "lon", member_dim = "member") } \arguments{ -\item{data}{A numerical array indicating the data to be used for the index -computation with latitude, longitude, start date, forecast month, and -member dimensions (in case of decadal predictions), with latitude, -longitude, year, month and member dimensions (in case of historical -simulations), or with latitude, longitude, year and month (in case of -observations or reanalyses). This data has to be provided, at least, over -the whole region needed to compute the index.} +\item{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 +provided, at least, over the whole region needed to compute the index.} \item{data_lats}{A numeric vector indicating the latitudes of the data.} @@ -73,9 +72,9 @@ dimension. The default value is 'member'. Only used if parameter 'type' is 'dcpp' or 'hist'.} } \value{ -An numerical array of the AMV index with the dimensions of: +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 the member (in case of historical simulations); or + 2) year and member (in case of historical simulations); or 3) year (in case of observations or reanalyses). } \description{ diff --git a/man/GMST.Rd b/man/GMST.Rd index 78508b5..3c9ace7 100644 --- a/man/GMST.Rd +++ b/man/GMST.Rd @@ -10,83 +10,122 @@ GMST(data_tas, data_tos, data_lats, data_lons, mask_sea_land, sea_value, type, year_dim = "year", month_dim = "month", member_dim = "member") } \arguments{ -\item{data_tas}{Surface air temperature data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), -with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). -This data has to be provided, at least, over the whole region needed to compute the index. The dimensions must be identical to those of data_tos.} - -\item{data_tos}{Sea surface temperature data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), -with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). -This data has to be provided, at least, over the whole region needed to compute the index. The dimensions must be identical to those of data_tas.} - -\item{data_lats}{An array with the latitudes of the data.} - -\item{data_lons}{An array with the longitudes of the data.} - -\item{mask_sea_land}{# mask for blending the surface air temperature over land and the sea surface temperature over the ocean. It must have lat_dim and lon_dim dimensions, and their lengths -must be equal to the length of data_lats and data_lons.} - -\item{sea_value}{Value of the sea grid points in the mask.} - -\item{type}{A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses).} - -\item{mask}{An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a -combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the -observational dataset for a fair comparison between the forecast system and the reference dataset.} - -\item{lat_dim}{A string with the name of the latitude dimension ('lat' by default).} - -\item{lon_dim}{A string with the name of the longitude dimension ('lon' by default).} - -\item{monini}{Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'.} - -\item{fmonth_dim}{A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'.} - -\item{sdate_dim}{A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'.} - -\item{indices_for_clim}{Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). -In case of type='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.} - -\item{year_dim}{A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'.} - -\item{month_dim}{A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'.} - -\item{member_dim}{A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'.} +\item{data_tas}{A numerical array indicating the surface air temperature data +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 provided, at least, +over the whole region needed to compute the index. The dimensions must be +identical to those of data_tos.} + +\item{data_tos}{A numerical array indicating the sea surface temperature data +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 provided, at least, +over the whole region needed to compute the index. The dimensions must be +identical to those of data_tas.} + +\item{data_lats}{A numeric vector indicating the latitudes of the data.} + +\item{data_lons}{A numeric vector indicating the longitudes of the data.} + +\item{mask_sea_land}{An array with dimensions [lat_dim = data_lats, lon_dim = +data_lons] for blending 'data_tas' and 'data_tos'.} + +\item{sea_value}{A numeric value indicating the sea grid points in +'mask_sea_land'.} + +\item{type}{A character string indicating the type of data ('dcpp' for +decadal predictions, 'hist' for historical simulations, or 'obs' for +observations or reanalyses).} + +\item{mask}{An array of a mask (with 0's in the grid points that have to be +masked) or NULL (i.e., no mask is used). This parameter allows to remove +the values over land in case the dataset is a combination of surface air +temperature over land and sea surface temperature over the ocean. Also, it +can be used to mask those grid points that are missing in the observational +dataset for a fair comparison between the forecast system and the reference +dataset. The default value is NULL.} + +\item{lat_dim}{A character string of the name of the latitude dimension. The +default value is 'lat'.} + +\item{lon_dim}{A character string of the name of the longitude dimension. The +default value is 'lon'.} + +\item{monini}{An integer indicating the month in which the forecast system is +initialized. Only used when parameter 'type' is 'dcpp'. The default value +is 11, i.e., initialized in November.} + +\item{fmonth_dim}{A character string indicating the name of the forecast +month dimension. Only used if parameter 'type' is 'dcpp'. The default value +is 'fmonth'.} + +\item{sdate_dim}{A character string indicating the name of the start date +dimension. Only used if parameter 'type' is 'dcpp'. The default value is +'sdate'.} + +\item{indices_for_clim}{A numeric vector of the indices of the years to +compute the climatology for calculating the anomalies, or NULL so the +climatology is calculated over the whole period. If the data are already +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.} + +\item{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'.} + +\item{month_dim}{A character string indicating the name of the month +dimension. The default value is 'month'. Only used if parameter 'type' is +'hist' or 'obs'.} + +\item{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'.} } \value{ -The GMST anomalies as function of the sdate, forecast year, and member (in case of decadal predictions); - as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +A numerical array of the GMST anomalies 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). } \description{ -The Global Mean Surface Temperature (GMST) anomalies are computed as the weighted-averaged -surface air temperature anomalies over land and sea surface temperature anomalies over the ocean. +The Global Mean Surface Temperature (GMST) anomalies are computed as the +weighted-averaged surface air temperature anomalies over land and sea surface +temperature anomalies over the ocean. } \examples{ ## Observations or reanalyses -obs_tas = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) -obs_tos = array(2:101, dim = c(year = 5, lat = 19, lon = 37, month = 12)) -mask_sea_land = array(c(1,0,1), dim = c(lat = 19, lon = 37)) -sea_value = 1 -lat = seq(-90, 90, 10) -lon = seq(0, 360, 10) -index_obs = GMST(data_tas = obs_tas, data_tos = obs_tos, data_lats = lat, data_lons = lon, type = 'obs', mask_sea_land = mask_sea_land, sea_value = sea_value) +obs_tas <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +obs_tos <- array(2:101, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +mask_sea_land <- array(c(1,0,1), dim = c(lat = 19, lon = 37)) +sea_value <- 1 +lat <- seq(-90, 90, 10) +lon <- seq(0, 360, 10) +index_obs <- GMST(data_tas = obs_tas, data_tos = obs_tos, data_lats = lat, data_lons = lon, type = 'obs', mask_sea_land = mask_sea_land, sea_value = sea_value) ## Historical simulations -hist_tas = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) -hist_tos = array(2:101, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) -mask_sea_land = array(c(1,0,1), dim = c(lat = 19, lon = 37)) -sea_value = 1 -lat = seq(-90, 90, 10) -lon = seq(0, 360, 10) -index_hist = GMST(data_tas = hist_tas, data_tos = hist_tos, data_lats = lat, data_lons = lon, type = 'hist', mask_sea_land = mask_sea_land, sea_value = sea_value) +hist_tas <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +hist_tos <- array(2:101, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +mask_sea_land <- array(c(1,0,1), dim = c(lat = 19, lon = 37)) +sea_value <- 1 +lat <- seq(-90, 90, 10) +lon <- seq(0, 360, 10) +index_hist <- GMST(data_tas = hist_tas, data_tos = hist_tos, data_lats = lat, data_lons = lon, type = 'hist', mask_sea_land = mask_sea_land, sea_value = sea_value) ## Decadal predictions -dcpp_tas = array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) -dcpp_tos = array(2:101, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) -mask_sea_land = array(c(1,0,1), dim = c(lat = 19, lon = 37)) -sea_value = 1 -lat = seq(-90, 90, 10) -lon = seq(0, 360, 10) -index_dcpp = GMST(data_tas = dcpp_tas, data_tos = dcpp_tos, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1, mask_sea_land = mask_sea_land, sea_value = sea_value) +dcpp_tas <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +dcpp_tos <- array(2:101, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +mask_sea_land <- array(c(1,0,1), dim = c(lat = 19, lon = 37)) +sea_value <- 1 +lat <- seq(-90, 90, 10) +lon <- seq(0, 360, 10) +index_dcpp <- GMST(data_tas = dcpp_tas, data_tos = dcpp_tos, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1, mask_sea_land = mask_sea_land, sea_value = sea_value) } \author{ diff --git a/man/GSAT.Rd b/man/GSAT.Rd index fda763e..bb43cfb 100644 --- a/man/GSAT.Rd +++ b/man/GSAT.Rd @@ -10,42 +10,72 @@ GSAT(data, data_lats, data_lons, type, lat_dim = "lat", lon_dim = "lon", member_dim = "member") } \arguments{ -\item{data}{Data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), -with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). -This data has to be provided, at least, over the whole region needed to compute the index.} - -\item{data_lats}{An array with the latitudes of the data.} - -\item{data_lons}{An array with the longitudes of the data.} - -\item{type}{A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses).} - -\item{lat_dim}{A string with the name of the latitude dimension ('lat' by default).} - -\item{lon_dim}{A string with the name of the longitude dimension ('lon' by default).} - -\item{mask}{An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a -combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the -observational dataset for a fair comparison between the forecast system and the reference dataset.} - -\item{monini}{Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'.} - -\item{fmonth_dim}{A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'.} - -\item{sdate_dim}{A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'.} - -\item{indices_for_clim}{Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). -In case of type='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.} - -\item{year_dim}{A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'.} - -\item{month_dim}{A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'.} - -\item{member_dim}{A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'.} +\item{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 +provided, at least, over the whole region needed to compute the index.} + +\item{data_lats}{A numeric vector indicating the latitudes of the data.} + +\item{data_lons}{A numeric vector indicating the longitudes of the data.} + +\item{type}{A character string indicating the type of data ('dcpp' for +decadal predictions, 'hist' for historical simulations, or 'obs' for +observations or reanalyses).} + +\item{lat_dim}{A character string of the name of the latitude dimension. The +default value is 'lat'.} + +\item{lon_dim}{A character string of the name of the longitude dimension. The +default value is 'lon'.} + +\item{mask}{An array of a mask (with 0's in the grid points that have to be +masked) or NULL (i.e., no mask is used). This parameter allows to remove +the values over land in case the dataset is a combination of surface air +temperature over land and sea surface temperature over the ocean. Also, it +can be used to mask those grid points that are missing in the observational +dataset for a fair comparison between the forecast system and the reference +dataset. The default value is NULL.} + +\item{monini}{An integer indicating the month in which the forecast system is +initialized. Only used when parameter 'type' is 'dcpp'. The default value +is 11, i.e., initialized in November.} + +\item{fmonth_dim}{A character string indicating the name of the forecast +month dimension. Only used if parameter 'type' is 'dcpp'. The default value +is 'fmonth'.} + +\item{sdate_dim}{A character string indicating the name of the start date +dimension. Only used if parameter 'type' is 'dcpp'. The default value is +'sdate'.} + +\item{indices_for_clim}{A numeric vector of the indices of the years to +compute the climatology for calculating the anomalies, or NULL so the +climatology is calculated over the whole period. If the data are already +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.} + +\item{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'.} + +\item{month_dim}{A character string indicating the name of the month +dimension. The default value is 'month'. Only used if parameter 'type' is +'hist' or 'obs'.} + +\item{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'.} } \value{ -The GSAT anomalies as function of the sdate, forecast year, and member (in case of decadal predictions); - as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +A numerical array of the GSAT anomalies 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). } \description{ The Global Surface Air Temperature (GSAT) anomalies are computed as the @@ -53,22 +83,22 @@ weighted-averaged surface air temperature anomalies over the global region. } \examples{ ## Observations or reanalyses -obs = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) -lat = seq(-90, 90, 10) -lon = seq(0, 360, 10) -index_obs = GSAT(data = obs, data_lats = lat, data_lons = lon, type = 'obs') +obs <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +lat <- seq(-90, 90, 10) +lon <- seq(0, 360, 10) +index_obs <- GSAT(data = obs, data_lats = lat, data_lons = lon, type = 'obs') ## Historical simulations -hist = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) -lat = seq(-90, 90, 10) -lon = seq(0, 360, 10) -index_hist = GSAT(data = hist, data_lats = lat, data_lons = lon, type = 'hist') +hist <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +lat <- seq(-90, 90, 10) +lon <- seq(0, 360, 10) +index_hist <- GSAT(data = hist, data_lats = lat, data_lons = lon, type = 'hist') ## Decadal predictions -dcpp = array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) -lat = seq(-90, 90, 10) -lon = seq(0, 360, 10) -index_dcpp = GSAT(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) +dcpp <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +lat <- seq(-90, 90, 10) +lon <- seq(0, 360, 10) +index_dcpp <- GSAT(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) } \author{ diff --git a/man/SPOD.Rd b/man/SPOD.Rd index 36ad341..09621c7 100644 --- a/man/SPOD.Rd +++ b/man/SPOD.Rd @@ -10,66 +10,98 @@ SPOD(data, data_lats, data_lons, type, lat_dim = "lat", lon_dim = "lon", member_dim = "member") } \arguments{ -\item{data}{Data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), -with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). -This data has to be provided, at least, over the whole region needed to compute the index.} - -\item{data_lats}{An array with the latitudes of the data.} - -\item{data_lons}{An array with the longitudes of the data.} - -\item{type}{A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses).} - -\item{lat_dim}{A string with the name of the latitude dimension ('lat' by default).} - -\item{lon_dim}{A string with the name of the longitude dimension ('lon' by default).} - -\item{mask}{An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a -combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the -observational dataset for a fair comparison between the forecast system and the reference dataset.} - -\item{monini}{Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'.} - -\item{fmonth_dim}{A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'.} - -\item{sdate_dim}{A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'.} - -\item{indices_for_clim}{Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). -In case of type='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.} - -\item{year_dim}{A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'.} - -\item{month_dim}{A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'.} - -\item{member_dim}{A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'.} +\item{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 +provided, at least, over the whole region needed to compute the index.} + +\item{data_lats}{A numeric vector indicating the latitudes of the data.} + +\item{data_lons}{A numeric vector indicating the longitudes of the data.} + +\item{type}{A character string indicating the type of data ('dcpp' for +decadal predictions, 'hist' for historical simulations, or 'obs' for +observations or reanalyses).} + +\item{lat_dim}{A character string of the name of the latitude dimension. The +default value is 'lat'.} + +\item{lon_dim}{A character string of the name of the longitude dimension. The +default value is 'lon'.} + +\item{mask}{An array of a mask (with 0's in the grid points that have to be +masked) or NULL (i.e., no mask is used). This parameter allows to remove +the values over land in case the dataset is a combination of surface air +temperature over land and sea surface temperature over the ocean. Also, it +can be used to mask those grid points that are missing in the observational +dataset for a fair comparison between the forecast system and the reference +dataset. The default value is NULL.} + +\item{monini}{An integer indicating the month in which the forecast system is +initialized. Only used when parameter 'type' is 'dcpp'. The default value +is 11, i.e., initialized in November.} + +\item{fmonth_dim}{A character string indicating the name of the forecast +month dimension. Only used if parameter 'type' is 'dcpp'. The default value +is 'fmonth'.} + +\item{sdate_dim}{A character string indicating the name of the start date +dimension. Only used if parameter 'type' is 'dcpp'. The default value is +'sdate'.} + +\item{indices_for_clim}{A numeric vector of the indices of the years to +compute the climatology for calculating the anomalies, or NULL so the +climatology is calculated over the whole period. If the data are already +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.} + +\item{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'.} + +\item{month_dim}{A character string indicating the name of the month +dimension. The default value is 'month'. Only used if parameter 'type' is +'hist' or 'obs'.} + +\item{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'.} } \value{ -The SPOD index as function of the sdate, forecast year, and member (in case of decadal predictions); - as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +A numerical array of the SPOD 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). } \description{ -The South Pacific Ocean Dipole (SPOD) index is related to the El Nino-Southern Oscillation (ENSO) and the Inderdecadal Pacific Oscillation (IPO). -The SPOD index is computed as the difference of weighted-averaged SST anomalies over 20ºS-48ºS, 165ºE-190ºE (NW pole) and the weighted-averaged SST anomalies over -44ºS-65ºS, 220ºE-260ºE (SE pole) (Saurral et al., 2020). +The South Pacific Ocean Dipole (SPOD) index is related to the El +Nino-Southern Oscillation (ENSO) and the Inderdecadal Pacific Oscillation +(IPO). The SPOD index is computed as the difference of weighted-averaged SST +anomalies over 20ºS-48ºS, 165ºE-190ºE (NW pole) and the weighted-averaged SST +anomalies over 44ºS-65ºS, 220ºE-260ºE (SE pole) (Saurral et al., 2020). } \examples{ ## Observations or reanalyses -obs = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) -lat = seq(-90, 90, 10) -lon = seq(0, 360, 10) -index_obs = SPOD(data = obs, data_lats = lat, data_lons = lon, type = 'obs') +obs <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) +lat <- seq(-90, 90, 10) +lon <- seq(0, 360, 10) +index_obs <- SPOD(data = obs, data_lats = lat, data_lons = lon, type = 'obs') ## Historical simulations -hist = array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) -lat = seq(-90, 90, 10) -lon = seq(0, 360, 10) -index_hist = SPOD(data = hist, data_lats = lat, data_lons = lon, type = 'hist') +hist <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) +lat <- seq(-90, 90, 10) +lon <- seq(0, 360, 10) +index_hist <- SPOD(data = hist, data_lats = lat, data_lons = lon, type = 'hist') ## Decadal predictions -dcpp = array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) -lat = seq(-90, 90, 10) -lon = seq(0, 360, 10) -index_dcpp = SPOD(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) +dcpp <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) +lat <- seq(-90, 90, 10) +lon <- seq(0, 360, 10) +index_dcpp <- SPOD(data = dcpp, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1) } \author{ diff --git a/man/TPI.Rd b/man/TPI.Rd index 0e5655b..d7938ff 100644 --- a/man/TPI.Rd +++ b/man/TPI.Rd @@ -10,46 +10,78 @@ TPI(data, data_lats, data_lons, index, type, lat_dim = "lat", month_dim = "month", member_dim = "member") } \arguments{ -\item{data}{Data to be used for the index computation with latitude, longitude, start date, forecast month, and member dimensions (in case of decadal predictions), -with latitude, longitude, year, month and member dimensions (in case of historical simulations), or with latitude, longitude, year and month (in case of observations or reanalyses). -This data has to be provided, at least, over the whole region needed to compute the index.} - -\item{data_lats}{An array with the latitudes of the data.} - -\item{data_lons}{An array with the longitudes of the data.} - -\item{type}{A string the the type of data ('dcpp' for decadal predictions, 'hist' for historical simulations, or 'obs' for observations or reanalyses).} - -\item{lat_dim}{A string with the name of the latitude dimension ('lat' by default).} - -\item{lon_dim}{A string with the name of the longitude dimension ('lon' by default).} - -\item{mask}{An array with a mask (with 0's in the grid points that have to be masked) or NULL (NULL by default, i.e., no mask is used). This parameter allows to remove the values over land in case the dataset is a -combination of surface air temperature over land and sea surface temperature over the ocean. Also, it can be used to mask those grid points that are missing in the -observational dataset for a fair comparison between the forecast system and the reference dataset.} - -\item{monini}{Month in which the forecast system is initialized (11 by default, i.e., initialized in November). Only used if type='dcpp'.} - -\item{fmonth_dim}{A string with the name of the forecast month dimension ('fmonth' by default). Only used if type='dcpp'.} - -\item{sdate_dim}{A string with the name of the start date dimension ('sdate' by default). Only used if type='dcpp'.} - -\item{indices_for_clim}{Indices of the years to compute the climatology. If NULL, the climatology is calculated over the whole period (NULL by default). -In case of type='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.} - -\item{year_dim}{A string with the name of the year dimension ('year' by default). Only used if type='hist' or type='obs'.} - -\item{month_dim}{A string with the name of the month dimension ('month' by default). Only used if type='hist' or type='obs'.} - -\item{member_dim}{A string with the name of the member dimension ('member' by default). Only used if type='dcpp' or type='hist'.} +\item{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 +provided, at least, over the whole region needed to compute the index.} + +\item{data_lats}{A numeric vector indicating the latitudes of the data.} + +\item{data_lons}{A numeric vector indicating the longitudes of the data.} + +\item{lat_dim}{A character string of the name of the latitude dimension. The +default value is 'lat'.} + +\item{lon_dim}{A character string of the name of the longitude dimension. The +default value is 'lon'.} + +\item{mask}{An array of a mask (with 0's in the grid points that have to be +masked) or NULL (i.e., no mask is used). This parameter allows to remove +the values over land in case the dataset is a combination of surface air +temperature over land and sea surface temperature over the ocean. Also, it +can be used to mask those grid points that are missing in the observational +dataset for a fair comparison between the forecast system and the reference +dataset. The default value is NULL.} + +\item{monini}{An integer indicating the month in which the forecast system is +initialized. Only used when parameter 'type' is 'dcpp'. The default value +is 11, i.e., initialized in November.} + +\item{fmonth_dim}{A character string indicating the name of the forecast +month dimension. Only used if parameter 'type' is 'dcpp'. The default value +is 'fmonth'.} + +\item{sdate_dim}{A character string indicating the name of the start date +dimension. Only used if parameter 'type' is 'dcpp'. The default value is +'sdate'.} + +\item{indices_for_clim}{A numeric vector of the indices of the years to +compute the climatology for calculating the anomalies, or NULL so the +climatology is calculated over the whole period. If the data are already +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.} + +\item{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'.} + +\item{month_dim}{A character string indicating the name of the month +dimension. The default value is 'month'. Only used if parameter 'type' is +'hist' or 'obs'.} + +\item{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'.} + +\item{A}{character string indicating the type of data ('dcpp' for +decadal predictions, 'hist' for historical simulations, or 'obs' for +observations or reanalyses).} } \value{ -The TPI as function of the sdate, forecast year, and member (in case of decadal predictions); - as function of the year and the member (in case of historical simulations); or as function of the year (in case of observations or reanalyses). +A numerical array of the TPI 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). } \description{ -The Tripole Index (TPI) for the Interdecadal Pacific Oscillation (IPO) is computed as the difference of weighted-averaged SST anomalies over -10ºS-10ºN, 170ºE-270ºE minus the mean of the weighted-averaged SST anomalies over 25ºN-45ºN, 140ºE-215ºE and 50ºS-15ºS, 150ºE-200ºE (Henley et al., 2015). +The Tripole Index (TPI) for the Interdecadal Pacific Oscillation (IPO) is +computed as the difference of weighted-averaged SST anomalies over 10ºS-10ºN, +170ºE-270ºE minus the mean of the weighted-averaged SST anomalies over +25ºN-45ºN, 140ºE-215ºE and 50ºS-15ºS, 150ºE-200ºE (Henley et al., 2015). } \examples{ ## Observations or reanalyses diff --git a/tests/testthat/test-GMST.R b/tests/testthat/test-GMST.R new file mode 100644 index 0000000..69d6556 --- /dev/null +++ b/tests/testthat/test-GMST.R @@ -0,0 +1,287 @@ +context("s2dv::GMST tests") + +############################################## + lat <- seq(-90, 90, 10) + lon <- seq(0, 360, 10) + set.seed(1) + mask <- array(sample(c(0, 1), 19*37, replace = T), dim = c(lat = 19, lon = 37)) + mask_sea_land <- array(c(1, 0, 1), dim = c(lat = 19, lon = 37)) + sea_value <- 1 + +## Observations or reanalyses + dat_tas1 <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) + dat_tos1 <- array(2:101, dim = c(year = 5, lat = 19, lon = 37, month = 12)) + + ## Historical simulations + dat_tas2 <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) + dat_tos2 <- array(2:101, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) + + ## Decadal predictions + dat_tas3 <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) + dat_tos3 <- array(2:101, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) + +############################################## +test_that("1. Input checks", { + + expect_error( + GMST(data_tas = c(), data_tos = c()), + "Parameter 'data_tas' and 'data_tos' cannot be NULL." + ) + expect_error( + GMST(data_tas = 'a', data_tos = dat_tos1), + "Parameter 'data_tas' and 'data_tos' must be a numeric array." + ) + expect_error( + GMST(data_tas = dat_tas1, data_tos = dat_tos2), + "The dimension of data_tas and data_tos must be identical." + ) + expect_error( + GMST(data_tas = dat_tas1, data_tos = dat_tos1, data_lats = c()), + "Parameter 'data_lats' must be a numeric vector." + ) + expect_error( + GMST(data_tas = dat_tas1, data_tos = dat_tos1, data_lats = lat, data_lons = c()), + "Parameter 'data_lons' must be a numeric vector." + ) + expect_error( + GMST(data_tas = dat_tas1, data_tos = dat_tos1, data_lats = lat, data_lons = lon, + mask_sea_land = TRUE), + "Parameter 'mask_sea_land' must be an array with dimensions \\[lat_dim, lon_dim\\]." + ) + expect_error( + GMST(data_tas = dat_tas1, data_tos = dat_tos1, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = TRUE), + "Parameter 'sea_value' must be a numeric value." + ) + expect_error( + GMST(data_tas = dat_tas1, data_tos = dat_tos1, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = TRUE), + "Parameter 'type' must be 'dcpp', 'hist', or 'obs'." + ) + expect_error( + GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp', + lat_dim = 2), + "Parameter 'lat_dim' must be a character string." + ) + expect_error( + GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp', + lat_dim = 'lats'), + "Parameter 'lat_dim' is not found in 'data_tas' or 'data_tos' dimension." + ) + expect_error( + GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp', + lon_dim = 2), + "Parameter 'lon_dim' must be a character string." + ) + expect_error( + GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp', + lon_dim = 'lons'), + "Parameter 'lon_dim' is not found in 'data_tas' or 'data_tos' dimension." + ) + expect_error( + GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = seq(-90, 90, 11), data_lons = lon, mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp'), + paste0("The latitude dimension of parameter 'data_tas' and 'data_tos'", + " must be the same length of parameter 'data_lats'.") + ) + expect_error( + GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = c(1:10), + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp'), + paste0("The longitude dimension of parameter 'data_tas' and 'data_tos'", " must be the same length of parameter 'data_lons'.") + ) + expect_error( + GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp', + mask = c(1, 0, 1, 0, 0)), + paste0("Parameter 'mask' must be NULL \\(no mask\\) or a numerical array ", + "with c\\(lat_dim, lon_dim\\) dimensions and 0 in those grid ", + "points that have to be masked.") + ) + expect_error( + GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp', + mask = array(rep(c(1, 0), 100), dim = c(lat = 20, lon = 5))), + paste0("Parameter 'mask' must be NULL \\(no mask\\) or a numerical array ", + "with c\\(lat_dim, lon_dim\\) dimensions and 0 in those grid ", + "points that have to be masked.") + ) + expect_error( + GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp', + monini = 0), + "Parameter 'monini' must be an integer from 1 to 12." + ) + expect_error( + GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp', + fmonth_dim = c('lat', 'lon')), + "Parameter 'fmonth_dim' must be a character string." + ) + expect_error( + GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp', + fmonth_dim = 'asd'), + "Parameter 'fmonth_dim' is not found in 'data_tas' or 'data_tos' dimension." + ) + expect_error( + GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp', + sdate_dim = 2), + "Parameter 'sdate_dim' must be a character string." + ) + expect_error( + GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp', + sdate_dim = 'date'), + "Parameter 'sdate_dim' is not found in 'data_tas' or 'data_tos' dimension." + ) + expect_error( + GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp', + indices_for_clim = TRUE), + paste0("The parameter 'indices_for_clim' must be a numeric vector ", + "or NULL to compute the anomalies based on the whole period, ", + "or FALSE if data are already anomalies") + ) + expect_error( + GMST(data_tas = dat_tas2, data_tos = dat_tos2, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'hist', + year_dim = 1), + "Parameter 'year_dim' must be a character string." + ) + expect_error( + GMST(data_tas = dat_tas2, data_tos = dat_tos2, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'hist', + year_dim = 'yr'), + "Parameter 'year_dim' is not found in 'data_tas' or 'data_tos' dimension." + ) + expect_error( + GMST(data_tas = dat_tas2, data_tos = dat_tos2, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'hist', + month_dim = 0), + "Parameter 'month_dim' must be a character string." + ) + expect_error( + GMST(data_tas = dat_tas2, data_tos = dat_tos2, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'hist', + month_dim = 'mth'), + "Parameter 'month_dim' is not found in 'data_tas' or 'data_tos' dimension." + ) + expect_error( + GMST(data_tas = dat_tas2, data_tos = dat_tos2, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'hist', + member_dim = 1), + "Parameter 'member_dim' must be a character string." + ) + expect_error( + GMST(data_tas = dat_tas2, data_tos = dat_tos2, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'hist', + member_dim = 'ens'), + "Parameter 'member_dim' is not found in 'data_tas' or 'data_tos' dimension." + ) + +}) + +############################################## +test_that("2. Output checks: dat1", { + res1_1 <- GMST(data_tas = dat_tas1, data_tos = dat_tos1, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'obs') + + expect_equal( + dim(res1_1), + c(year = 5) + ) + expect_equal( + res1_1, + array(seq(-2, 2, 1), dim = c(year = 5)) + ) + + res1_2 <- GMST(data_tas = dat_tas1, data_tos = dat_tos1, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'obs', + indices_for_clim = F) + expect_equal( + mean(res1_2), + 51.09007, + tolerance = 0.00001 + ) + +}) +############################################## +test_that("3. Output checks: dat2", { + res2_1 <- GMST(data_tas = dat_tas2, data_tos = dat_tos2, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'hist') + expect_equal( + dim(res2_1), + c(year = 5, member = 5) + ) + expect_equal( + mean(res2_1), + 0 + ) + expect_equal( + range(res2_1), + c(-2.028458, 2.028458), + tolerance = 0.00001 + ) + + res2_2 <- GMST(data_tas = dat_tas2, data_tos = dat_tos2, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'hist', + indices_for_clim = 2:3) + expect_equal( + mean(res2_2), + 0.5 + ) + + res2_3 <- GMST(data_tas = dat_tas2, data_tos = dat_tos2, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'hist', + mask = mask) + expect_equal( + mean(res2_3)*10^16, + 2.842084, + tolerance = 0.00001 + ) + +}) +############################################## +test_that("4. Output checks: dat3", { + res3_1 <- GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp') + expect_equal( + dim(res3_1), + c(fyear = 1, sdate = 5, member = 5) + ) + expect_equal( + mean(res3_1)*10^16, + -2.842084, + tolerance = 0.00001 + ) + + res3_2 <- GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp', + indices_for_clim = 2:3) + expect_equal( + mean(res3_2), + 0.5 + ) + + res3_3 <- GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp', + monini = 10) + expect_equal( + mean(res3_3), + 0 + ) + + res3_4 <- GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, + mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp', + mask = mask) + expect_equal( + mean(res3_4)*10^16, + -8.52686, + tolerance = 0.00001 + ) + +}) diff --git a/tests/testthat/test-GSAT.R b/tests/testthat/test-GSAT.R new file mode 100644 index 0000000..bd53e8a --- /dev/null +++ b/tests/testthat/test-GSAT.R @@ -0,0 +1,240 @@ +context("s2dv::GSAT tests") + +############################################## + lat <- seq(-90, 90, 10) + lon <- seq(0, 360, 10) + set.seed(1) + mask <- array(sample(c(0, 1), 19*37, replace = T), dim = c(lat = 19, lon = 37)) + +## Observations or reanalyses + dat1 <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) + + ## Historical simulations + dat2 <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) + + ## Decadal predictions + dat3 <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) + +############################################## +test_that("1. Input checks", { + + expect_error( + GSAT(data = c()), + "Parameter 'data' cannot be NULL." + ) + expect_error( + GSAT(data = 'a'), + "Parameter 'data' must be a numeric array." + ) + expect_error( + GSAT(data = dat1, data_lats = c()), + "Parameter 'data_lats' must be a numeric vector." + ) + expect_error( + GSAT(data = dat1, data_lats = lat, data_lons = c()), + "Parameter 'data_lons' must be a numeric vector." + ) + expect_error( + GSAT(data = dat1, data_lats = lat, data_lons = lon, type = TRUE), + "Parameter 'type' must be 'dcpp', 'hist', or 'obs'." + ) + expect_error( + GSAT(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + lat_dim = 2), + "Parameter 'lat_dim' must be a character string." + ) + expect_error( + GSAT(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + lat_dim = 'lats'), + "Parameter 'lat_dim' is not found in 'data' dimension." + ) + expect_error( + GSAT(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + lon_dim = 2), + "Parameter 'lon_dim' must be a character string." + ) + expect_error( + GSAT(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + lon_dim = 'lons'), + "Parameter 'lon_dim' is not found in 'data' dimension." + ) + expect_error( + GSAT(data = dat3, data_lats = seq(-90, 90, 11), data_lons = lon, type = 'dcpp'), + paste0("The latitude dimension of parameter 'data' must be the same", + " length of parameter 'data_lats'.") + ) + expect_error( + GSAT(data = dat3, data_lats = lat, data_lons = c(1:10), type = 'dcpp'), + paste0("The longitude dimension of parameter 'data' must be the same", + " length of parameter 'data_lons'.") + ) + expect_error( + GSAT(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + mask = c(1, 0, 1, 0, 0)), + paste0("Parameter 'mask' must be NULL \\(no mask\\) or a numerical array ", + "with c\\(lat_dim, lon_dim\\) dimensions and 0 in those grid ", + "points that have to be masked.") + ) + expect_error( + GSAT(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + mask = array(rep(c(1, 0), 100), dim = c(lat = 20, lon = 5))), + paste0("Parameter 'mask' must be NULL \\(no mask\\) or a numerical array ", + "with c\\(lat_dim, lon_dim\\) dimensions and 0 in those grid ", + "points that have to be masked.") + ) + expect_error( + GSAT(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + monini = 0), + "Parameter 'monini' must be an integer from 1 to 12." + ) + expect_error( + GSAT(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + fmonth_dim = c('lat', 'lon')), + "Parameter 'fmonth_dim' must be a character string." + ) + expect_error( + GSAT(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + fmonth_dim = 'asd'), + "Parameter 'fmonth_dim' is not found in 'data' dimension." + ) + expect_error( + GSAT(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + sdate_dim = 2), + "Parameter 'sdate_dim' must be a character string." + ) + expect_error( + GSAT(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + sdate_dim = 'date'), + "Parameter 'sdate_dim' is not found in 'data' dimension." + ) + expect_error( + GSAT(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + indices_for_clim = TRUE), + paste0("The parameter 'indices_for_clim' must be a numeric vector ", + "or NULL to compute the anomalies based on the whole period, ", + "or FALSE if data are already anomalies") + ) + expect_error( + GSAT(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + year_dim = 1), + "Parameter 'year_dim' must be a character string." + ) + expect_error( + GSAT(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + year_dim = 'yr'), + "Parameter 'year_dim' is not found in 'data' dimension." + ) + expect_error( + GSAT(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + month_dim = 0), + "Parameter 'month_dim' must be a character string." + ) + expect_error( + GSAT(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + month_dim = 'mth'), + "Parameter 'month_dim' is not found in 'data' dimension." + ) + expect_error( + GSAT(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + member_dim = 1), + "Parameter 'member_dim' must be a character string." + ) + expect_error( + GSAT(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + member_dim = 'ens'), + "Parameter 'member_dim' is not found in 'data' dimension." + ) + +}) + +############################################## +test_that("2. Output checks: dat1", { + res1_1 <- GSAT(data = dat1, data_lats = lat, data_lons = lon, type = 'obs') + + expect_equal( + dim(res1_1), + c(year = 5) + ) + expect_equal( + mean(res1_1)*10^15, + 1.421042, + tolerance = 0.00001 + ) + + res1_2 <- GSAT(data = dat1, data_lats = lat, data_lons = lon, type = 'obs', + indices_for_clim = F) + expect_equal( + mean(res1_2), + 50.42349, + tolerance = 0.00001 + ) + +}) +############################################## +test_that("3. Output checks: dat2", { + res2_1 <- GSAT(data = dat2, data_lats = lat, data_lons = lon, type = 'hist') + expect_equal( + dim(res2_1), + c(year = 5, member = 5) + ) + expect_equal( + mean(res2_1)*10^16, + 2.841867, + tolerance = 0.00001 + ) + + res2_2 <- GSAT(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + indices_for_clim = 2:3) + expect_equal( + mean(res2_2), + 0.5, + tolerance = 0.00001 + ) + + res2_3 <- GSAT(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + mask = mask) + expect_equal( + mean(res2_3)*10^15, + 1.136842, + tolerance = 0.00001 + ) + +}) +############################################## +test_that("4. Output checks: dat3", { + res3_1 <- GSAT(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp') + expect_equal( + dim(res3_1), + c(fyear = 1, sdate = 5, member = 5) + ) + expect_equal( + mean(res3_1)*10^16, + -5.684776, + tolerance = 0.00001 + ) + + res3_2 <- GSAT(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + indices_for_clim = 2:3) + expect_equal( + mean(res3_2), + 0.5, + tolerance = 0.00001 + ) + + res3_3 <- GSAT(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + monini = 5) + expect_equal( + mean(res3_3)*10^16, + 2.841867, + tolerance = 0.00001 + ) + + res3_4 <- GSAT(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + mask = mask) + expect_equal( + mean(res3_4), + 0, + tolerance = 0.00001 + ) + +}) diff --git a/tests/testthat/test-SPOD.R b/tests/testthat/test-SPOD.R new file mode 100644 index 0000000..6d23a59 --- /dev/null +++ b/tests/testthat/test-SPOD.R @@ -0,0 +1,240 @@ +context("s2dv::SPOD tests") + +############################################## + lat <- seq(-90, 90, 10) + lon <- seq(0, 360, 10) + set.seed(1) + mask <- array(sample(c(0, 1), 19*37, replace = T), dim = c(lat = 19, lon = 37)) + +## Observations or reanalyses + dat1 <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) + + ## Historical simulations + dat2 <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) + + ## Decadal predictions + dat3 <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) + +############################################## +test_that("1. Input checks", { + + expect_error( + SPOD(data = c()), + "Parameter 'data' cannot be NULL." + ) + expect_error( + SPOD(data = 'a'), + "Parameter 'data' must be a numeric array." + ) + expect_error( + SPOD(data = dat1, data_lats = c()), + "Parameter 'data_lats' must be a numeric vector." + ) + expect_error( + SPOD(data = dat1, data_lats = lat, data_lons = c()), + "Parameter 'data_lons' must be a numeric vector." + ) + expect_error( + SPOD(data = dat1, data_lats = lat, data_lons = lon, type = TRUE), + "Parameter 'type' must be 'dcpp', 'hist', or 'obs'." + ) + expect_error( + SPOD(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + lat_dim = 2), + "Parameter 'lat_dim' must be a character string." + ) + expect_error( + SPOD(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + lat_dim = 'lats'), + "Parameter 'lat_dim' is not found in 'data' dimension." + ) + expect_error( + SPOD(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + lon_dim = 2), + "Parameter 'lon_dim' must be a character string." + ) + expect_error( + SPOD(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + lon_dim = 'lons'), + "Parameter 'lon_dim' is not found in 'data' dimension." + ) + expect_error( + SPOD(data = dat3, data_lats = seq(-90, 90, 11), data_lons = lon, type = 'dcpp'), + paste0("The latitude dimension of parameter 'data' must be the same", + " length of parameter 'data_lats'.") + ) + expect_error( + SPOD(data = dat3, data_lats = lat, data_lons = c(1:10), type = 'dcpp'), + paste0("The longitude dimension of parameter 'data' must be the same", + " length of parameter 'data_lons'.") + ) + expect_error( + SPOD(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + mask = c(1, 0, 1, 0, 0)), + paste0("Parameter 'mask' must be NULL \\(no mask\\) or a numerical array ", + "with c\\(lat_dim, lon_dim\\) dimensions and 0 in those grid ", + "points that have to be masked.") + ) + expect_error( + SPOD(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + mask = array(rep(c(1, 0), 100), dim = c(lat = 20, lon = 5))), + paste0("Parameter 'mask' must be NULL \\(no mask\\) or a numerical array ", + "with c\\(lat_dim, lon_dim\\) dimensions and 0 in those grid ", + "points that have to be masked.") + ) + expect_error( + SPOD(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + monini = 0), + "Parameter 'monini' must be an integer from 1 to 12." + ) + expect_error( + SPOD(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + fmonth_dim = c('lat', 'lon')), + "Parameter 'fmonth_dim' must be a character string." + ) + expect_error( + SPOD(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + fmonth_dim = 'asd'), + "Parameter 'fmonth_dim' is not found in 'data' dimension." + ) + expect_error( + SPOD(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + sdate_dim = 2), + "Parameter 'sdate_dim' must be a character string." + ) + expect_error( + SPOD(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + sdate_dim = 'date'), + "Parameter 'sdate_dim' is not found in 'data' dimension." + ) + expect_error( + SPOD(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + indices_for_clim = TRUE), + paste0("The parameter 'indices_for_clim' must be a numeric vector ", + "or NULL to compute the anomalies based on the whole period, ", + "or FALSE if data are already anomalies") + ) + expect_error( + SPOD(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + year_dim = 1), + "Parameter 'year_dim' must be a character string." + ) + expect_error( + SPOD(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + year_dim = 'yr'), + "Parameter 'year_dim' is not found in 'data' dimension." + ) + expect_error( + SPOD(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + month_dim = 0), + "Parameter 'month_dim' must be a character string." + ) + expect_error( + SPOD(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + month_dim = 'mth'), + "Parameter 'month_dim' is not found in 'data' dimension." + ) + expect_error( + SPOD(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + member_dim = 1), + "Parameter 'member_dim' must be a character string." + ) + expect_error( + SPOD(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + member_dim = 'ens'), + "Parameter 'member_dim' is not found in 'data' dimension." + ) + +}) + +############################################## +test_that("2. Output checks: dat1", { + res1_1 <- SPOD(data = dat1, data_lats = lat, data_lons = lon, type = 'obs') + + expect_equal( + dim(res1_1), + c(year = 5) + ) + expect_equal( + mean(res1_1)*10^16, + 1.776357, + tolerance = 0.00001 + ) + + res1_2 <- SPOD(data = dat1, data_lats = lat, data_lons = lon, type = 'obs', + indices_for_clim = F) + expect_equal( + mean(res1_2), + 6.073599, + tolerance = 0.00001 + ) + +}) +############################################## +test_that("3. Output checks: dat2", { + res2_1 <- SPOD(data = dat2, data_lats = lat, data_lons = lon, type = 'hist') + expect_equal( + dim(res2_1), + c(year = 5, member = 5) + ) + expect_equal( + mean(res2_1)*10^16, + 2.754588, + tolerance = 0.00001 + ) + + res2_2 <- SPOD(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + indices_for_clim = 2:3) + expect_equal( + mean(res2_2)*10^16, + 2.645063, + tolerance = 0.00001 + ) + + res2_3 <- SPOD(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + mask = mask) + expect_equal( + mean(res2_3)*10^17, + 5.3553, + tolerance = 0.00001 + ) + +}) +############################################## +test_that("4. Output checks: dat3", { + res3_1 <- SPOD(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp') + expect_equal( + dim(res3_1), + c(fyear = 1, sdate = 5, member = 5) + ) + expect_equal( + mean(res3_1)*10^17, + -3.552714, + tolerance = 0.00001 + ) + + res3_2 <- SPOD(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + indices_for_clim = 2:3) + expect_equal( + mean(res3_2)*10^16, + 1.421085, + tolerance = 0.00001 + ) + + res3_3 <- SPOD(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + monini = 5) + expect_equal( + mean(res3_3)*10^16, + 2.087219, + tolerance = 0.00001 + ) + + res3_4 <- SPOD(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + mask = mask) + expect_equal( + mean(res3_4)*10^17, + 7.549517, + tolerance = 0.00001 + ) + +}) diff --git a/tests/testthat/test-TPI.R b/tests/testthat/test-TPI.R new file mode 100644 index 0000000..3b584c6 --- /dev/null +++ b/tests/testthat/test-TPI.R @@ -0,0 +1,245 @@ +context("s2dv::TPI tests") + +############################################## + lat <- seq(-90, 90, 10) + lon <- seq(0, 360, 10) + set.seed(1) + mask <- array(sample(c(0, 1), 19*37, replace = T), dim = c(lat = 19, lon = 37)) + +## Observations or reanalyses + dat1 <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12)) + + ## Historical simulations + dat2 <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5)) + + ## Decadal predictions + dat3 <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5)) + +############################################## +test_that("1. Input checks", { + + expect_error( + TPI(data = c()), + "Parameter 'data' cannot be NULL." + ) + expect_error( + TPI(data = 'a'), + "Parameter 'data' must be a numeric array." + ) + expect_error( + TPI(data = dat1, data_lats = c()), + "Parameter 'data_lats' must be a numeric vector." + ) + expect_error( + TPI(data = dat1, data_lats = lat, data_lons = c()), + "Parameter 'data_lons' must be a numeric vector." + ) + expect_error( + TPI(data = dat1, data_lats = lat, data_lons = lon, type = TRUE), + "Parameter 'type' must be 'dcpp', 'hist', or 'obs'." + ) + expect_error( + TPI(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + lat_dim = 2), + "Parameter 'lat_dim' must be a character string." + ) + expect_error( + TPI(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + lat_dim = 'lats'), + "Parameter 'lat_dim' is not found in 'data' dimension." + ) + expect_error( + TPI(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + lon_dim = 2), + "Parameter 'lon_dim' must be a character string." + ) + expect_error( + TPI(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + lon_dim = 'lons'), + "Parameter 'lon_dim' is not found in 'data' dimension." + ) + expect_error( + TPI(data = dat3, data_lats = seq(-90, 90, 11), data_lons = lon, type = 'dcpp'), + paste0("The latitude dimension of parameter 'data' must be the same", + " length of parameter 'data_lats'.") + ) + expect_error( + TPI(data = dat3, data_lats = lat, data_lons = c(1:10), type = 'dcpp'), + paste0("The longitude dimension of parameter 'data' must be the same", + " length of parameter 'data_lons'.") + ) + expect_error( + TPI(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + mask = c(1, 0, 1, 0, 0)), + paste0("Parameter 'mask' must be NULL \\(no mask\\) or a numerical array ", + "with c\\(lat_dim, lon_dim\\) dimensions and 0 in those grid ", + "points that have to be masked.") + ) + expect_error( + TPI(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + mask = array(rep(c(1, 0), 100), dim = c(lat = 20, lon = 5))), + paste0("Parameter 'mask' must be NULL \\(no mask\\) or a numerical array ", + "with c\\(lat_dim, lon_dim\\) dimensions and 0 in those grid ", + "points that have to be masked.") + ) + expect_error( + TPI(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + monini = 0), + "Parameter 'monini' must be an integer from 1 to 12." + ) + expect_error( + TPI(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + fmonth_dim = c('lat', 'lon')), + "Parameter 'fmonth_dim' must be a character string." + ) + expect_error( + TPI(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + fmonth_dim = 'asd'), + "Parameter 'fmonth_dim' is not found in 'data' dimension." + ) + expect_error( + TPI(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + sdate_dim = 2), + "Parameter 'sdate_dim' must be a character string." + ) + expect_error( + TPI(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + sdate_dim = 'date'), + "Parameter 'sdate_dim' is not found in 'data' dimension." + ) + expect_error( + TPI(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + indices_for_clim = TRUE), + paste0("The parameter 'indices_for_clim' must be a numeric vector ", + "or NULL to compute the anomalies based on the whole period, ", + "or FALSE if data are already anomalies") + ) + expect_error( + TPI(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + year_dim = 1), + "Parameter 'year_dim' must be a character string." + ) + expect_error( + TPI(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + year_dim = 'yr'), + "Parameter 'year_dim' is not found in 'data' dimension." + ) + expect_error( + TPI(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + month_dim = 0), + "Parameter 'month_dim' must be a character string." + ) + expect_error( + TPI(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + month_dim = 'mth'), + "Parameter 'month_dim' is not found in 'data' dimension." + ) + expect_error( + TPI(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + member_dim = 1), + "Parameter 'member_dim' must be a character string." + ) + expect_error( + TPI(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + member_dim = 'ens'), + "Parameter 'member_dim' is not found in 'data' dimension." + ) + +}) + +############################################## +test_that("2. Output checks: dat1", { + res1_1 <- TPI(data = dat1, data_lats = lat, data_lons = lon, type = 'obs') + + expect_equal( + dim(res1_1), + c(year = 5) + ) + expect_equal( + mean(res1_1), + 0, + tolerance = 0.00001 + ) + expect_equal( + res1_1*10^16, + array(c(0.000000, -12.21245, -5.551115, 12.21245, 5.551115), dim = c(year = 5)), + tolerance = 0.00001 + ) + + res1_2 <- TPI(data = dat1, data_lats = lat, data_lons = lon, type = 'obs', + indices_for_clim = F) + expect_equal( + mean(res1_2), + 0.9376566, + tolerance = 0.00001 + ) + +}) +############################################## +test_that("3. Output checks: dat2", { + res2_1 <- TPI(data = dat2, data_lats = lat, data_lons = lon, type = 'hist') + expect_equal( + dim(res2_1), + c(year = 5, member = 5) + ) + expect_equal( + mean(res2_1)*10^17, + -2.220099, + tolerance = 0.00001 + ) + + res2_2 <- TPI(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + indices_for_clim = 2:3) + expect_equal( + mean(res2_2)*10^16, + 3.840678, + tolerance = 0.00001 + ) + + res2_3 <- TPI(data = dat2, data_lats = lat, data_lons = lon, type = 'hist', + mask = mask) + expect_equal( + mean(res2_3)*10^17, + -2.288135, + tolerance = 0.00001 + ) + +}) +############################################## +test_that("4. Output checks: dat3", { + res3_1 <- TPI(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp') + expect_equal( + dim(res3_1), + c(fyear = 1, sdate = 5, member = 5) + ) + expect_equal( + mean(res3_1)*10^17, + 1.776357, + tolerance = 0.00001 + ) + + res3_2 <- TPI(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + indices_for_clim = 2:3) + expect_equal( + mean(res3_2)*10^16, + 3.064216, + tolerance = 0.00001 + ) + + res3_3 <- TPI(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + monini = 5) + expect_equal( + mean(res3_3)*10^17, + -2.664535, + tolerance = 0.00001 + ) + + res3_4 <- TPI(data = dat3, data_lats = lat, data_lons = lon, type = 'dcpp', + mask = mask) + expect_equal( + mean(res3_4)*10^17, + -1.360023, + tolerance = 0.00001 + ) + +}) -- GitLab From 9e70f6a1471b7a4bd86ba6dadcb8e9b4c4fc9695 Mon Sep 17 00:00:00 2001 From: aho Date: Thu, 5 Nov 2020 10:16:43 +0100 Subject: [PATCH 11/11] Adjust the input check order --- R/GMST.R | 59 +++++++++++++++++++------------------- tests/testthat/test-GMST.R | 5 ++-- 2 files changed, 33 insertions(+), 31 deletions(-) diff --git a/R/GMST.R b/R/GMST.R index c827951..5e6f4af 100644 --- a/R/GMST.R +++ b/R/GMST.R @@ -154,8 +154,12 @@ GMST <- function(data_tas, data_tos, data_lats, data_lons, mask_sea_land, sea_va stop(paste0("The longitude dimension of parameter 'data_tas' and 'data_tos'", " must be the same length of parameter 'data_lons'.")) } + # sea_value + if (!is.numeric(sea_value) | length(sea_value) != 1) { + stop("Parameter 'sea_value' must be a numeric value.") + } # mask_sea_land - if (!is.array(mask_sea_land)){ + if (!is.array(mask_sea_land)) { stop("Parameter 'mask_sea_land' must be an array with dimensions [lat_dim, lon_dim].") } else if (!identical(names(dim(mask_sea_land)), c(lat_dim, lon_dim))) { stop("Parameter 'mask_sea_land' must be an array with dimensions [lat_dim, lon_dim].") @@ -163,40 +167,14 @@ GMST <- function(data_tas, data_tos, data_lats, data_lons, mask_sea_land, sea_va c(length(data_lats), length(data_lons)))) { stop("Parameter 'mask_sea_land' dimensions must be equal to the length of 'data_lats' and 'data_lons'.") } - ## combination of tas and tos (data) - mask_tas_tos <- function(data_tas, data_tos, mask_sea_land, sea_value) { - data <- data_tas - data[mask_sea_land == sea_value] <- data_tos[mask_sea_land == sea_value] - return(data) - } - mask_sea_land <- s2dv::Reorder(data = mask_sea_land, order = c(lat_dim,lon_dim)) - data <- multiApply::Apply(data = list(data_tas, data_tos), - target_dims = c(lat_dim, lon_dim), - fun = mask_tas_tos, mask_sea_land = mask_sea_land, - sea_value = sea_value)$output1 - data <- drop(data) - - # sea_value - if (!is.numeric(sea_value) | length(sea_value) != 1) { - stop("Parameter 'sea_value' must be a numeric value.") - } # type if (!type %in% c('dcpp', 'hist', 'obs')) { stop("Parameter 'type' must be 'dcpp', 'hist', or 'obs'.") } # mask if (!is.null(mask)) { - if (is.array(mask) & identical(names(dim(mask)), c(lat_dim,lon_dim)) & - identical(as.integer(dim(mask)), c(length(data_lats), length(data_lons)))) { - ## To mask those grid point that are missing in the observations - mask <- s2dv::Reorder(data = mask, order = c(lat_dim, lon_dim)) - fun_mask <- function(data, mask) { - data[mask == 0] <- NA - return(data) - } - data <- multiApply::Apply(data = data, target_dims = c(lat_dim, lon_dim), - fun = fun_mask, mask = mask)$output1 - } else { + if (!is.array(mask) | !identical(names(dim(mask)), c(lat_dim,lon_dim)) | + !identical(as.integer(dim(mask)), c(length(data_lats), length(data_lons)))) { stop(paste0("Parameter 'mask' must be NULL (no mask) or a numerical array ", "with c(lat_dim, lon_dim) dimensions and 0 in those grid ", "points that have to be masked.")) @@ -264,8 +242,31 @@ GMST <- function(data_tas, data_tos, data_lats, data_lons, mask_sea_land, sea_va } } + ## combination of tas and tos (data) + mask_tas_tos <- function(data_tas, data_tos, mask_sea_land, sea_value) { + data <- data_tas + data[mask_sea_land == sea_value] <- data_tos[mask_sea_land == sea_value] + return(data) + } + mask_sea_land <- s2dv::Reorder(data = mask_sea_land, order = c(lat_dim,lon_dim)) + data <- multiApply::Apply(data = list(data_tas, data_tos), + target_dims = c(lat_dim, lon_dim), + fun = mask_tas_tos, mask_sea_land = mask_sea_land, + sea_value = sea_value)$output1 + data <- drop(data) rm(data_tas, data_tos) + ## To mask those grid point that are missing in the observations + if (!is.null(mask)) { + mask <- s2dv::Reorder(data = mask, order = c(lat_dim, lon_dim)) + fun_mask <- function(data, mask) { + data[mask == 0] <- NA + return(data) + } + data <- multiApply::Apply(data = data, target_dims = c(lat_dim, lon_dim), + fun = fun_mask, mask = mask)$output1 + } + data <- ClimProjDiags::WeightedMean(data = data, lon = data_lons, lat = data_lats, region = NULL, londim = which(names(dim(data)) == lon_dim), diff --git a/tests/testthat/test-GMST.R b/tests/testthat/test-GMST.R index 69d6556..d6e3a4f 100644 --- a/tests/testthat/test-GMST.R +++ b/tests/testthat/test-GMST.R @@ -45,7 +45,7 @@ test_that("1. Input checks", { ) expect_error( GMST(data_tas = dat_tas1, data_tos = dat_tos1, data_lats = lat, data_lons = lon, - mask_sea_land = TRUE), + sea_value = sea_value, mask_sea_land = TRUE), "Parameter 'mask_sea_land' must be an array with dimensions \\[lat_dim, lon_dim\\]." ) expect_error( @@ -90,7 +90,8 @@ test_that("1. Input checks", { expect_error( GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = c(1:10), mask_sea_land = mask_sea_land, sea_value = sea_value, type = 'dcpp'), - paste0("The longitude dimension of parameter 'data_tas' and 'data_tos'", " must be the same length of parameter 'data_lons'.") + paste0("The longitude dimension of parameter 'data_tas' and 'data_tos'", + " must be the same length of parameter 'data_lons'.") ) expect_error( GMST(data_tas = dat_tas3, data_tos = dat_tos3, data_lats = lat, data_lons = lon, -- GitLab