Skip to content
GitLab
Commits (12)
Hide whitespace changes
Inline Side-by-side
DESCRIPTION
View file @ dca0f085
Package: ClimProjDiags
Title: Set of Tools to Compute Various Climate Indices
Version: 0.1.1
Version: 0.1.3
Authors@R: c(
person("BSC-CNS", role = c("aut", "cph")),
person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-8568-3071")),
person("An-Chi", "Ho", , "an.ho@bsc.es", role = "ctb"),
person("Nicolau", "Manubens", , "nicolau.manubens@bsc.es", role = "ctb"),
person("Alasdair", "Hunter", , "alasdair.hunter@bsc.es", role = "aut"),
person("Louis-Philippe", "Caron", , "louis-philippe.caron@bsc.es", role = "ctb"))
......@@ -26,7 +27,6 @@ License: Apache License 2.0
URL: https://earth.bsc.es/gitlab/es/ClimProjDiags
BugReports: https://earth.bsc.es/gitlab/es/ClimProjDiags/-/issues
Encoding: UTF-8
LazyData: true
RoxygenNote: 5.0.0
Suggests:
knitr,
......
NAMESPACE
View file @ dca0f085
# Generated by roxygen2: do not edit by hand
export(AnoAgree)
export(ArrayToList)
export(Climdex)
export(CombineIndices)
export(DTRIndicator)
......
R/ArrayToList.R 0 → 100644
View file @ dca0f085
#' Split an array into list by a given array dimension
#'
#'@description This function splits an array into a list as required by PlotLayout function from package "s2dv" when parameter 'special_args' is used. The function ArrayToList allows to add names to the elements of the list in two different levels, the 'list' or the 'sublist'.
#'
#'@param data A multidimensional array.
#'@param dim A character string indicating the name of the dimension to split or an integer indicating the position of the dimension.
#'@param level A string character 'list' or 'sublist' indicating if it should be a list or a sublist. By default it creates a list.
#'@param names A vector of character strings to name the list (if it is a single string, it would be reused) or a single character string to name the elements in the sublist.
#'
#'@return A list of arrays of the length of the dimension set in parameter 'dim'.
#'
#'@examples
#'data <- array(1:240, c(month = 12, member = 5, time = 4))
#'# Create a list:
#'datalist <- ArrayToList(data, dim = 'month', level = 'list', names = month.name)
#'class(datalist)
#'class(datalist[[1]])
#'str(datalist)
#'# Create a sublist:
#'datalist <- ArrayToList(data, dim = 'month', level = 'sublist', names = 'dots')
#'class(datalist)
#'class(datalist[[1]])
#'class(datalist[[1]][[1]])
#'str(datalist)
#'@seealso \link[s2dv]{PlotLayout}
#'@export
ArrayToList <- function(data, dim, level = 'list', names = NULL) {
if (is.null(dim(data))) {
stop("Parameter 'data' must be an array or matrix.")
}
if (length(dim) > 1) {
warning("The lenght of parameter 'dim' is greater than 1 and only the",
" first element is used.")
dim <- dim[1]
}
if (is.numeric(dim) & dim > length(dim(data))) {
stop("Parameter 'dim' cannot exceed the length of 'data' dimensions.")
}
if (is.character(dim)) {
dim <- which(names(dim(data)) == dim)
if (identical(dim, integer(0))) {
stop("Parameter 'dim' is not found in the dimension names.")
}
}
if( !is.character(names)) {
stop("Parameter 'names' is not class character.")
}
datalist <- asplit(data, dim)
attributes(datalist) <- NULL
if (level == 'list') {
if (length(names) == 1) {
names <- rep(names, dim(data)[dim])
} else {
if (length(names) != dim(data)[dim]) {
stop("The length of parameter names could be 1 or equal to the dimension to split.")
}
}
names(datalist) <- names
} else if (level == 'sublist') {
if (length(names) > 1) {
warning("Parameter 'names' should be a single character string to name the sublist.",
"Only the first element is used.")
names <- names[1]
}
datalist <- lapply(datalist, function(x) {
res <- list(x)
names(res) <- names
return(res)})
} else {
stop("Parameter 'level' should be 'list' nor 'sublist'.")
}
return(datalist)
}
R/WeightedMean.R
View file @ dca0f085
......@@ -103,9 +103,13 @@ WeightedMean <- function(data, lon, lat, region = NULL, mask = NULL, londim = NU
mask <- aux$mask
}
}
wtmean <- aaply(data, .margins = dims[c(-londim, -latdim)],
.fun = .WeightedMean, lon, lat, mask, .drop = FALSE)
dim(wtmean) <- dim(data)[-c(latdim,londim)]
if (length(dim(data)) > 2) {
wtmean <- aaply(data, .margins = dims[c(-londim, -latdim)],
.fun = .WeightedMean, lon, lat, mask, .drop = FALSE)
dim(wtmean) <- dim(data)[-c(latdim,londim)]
} else {
wtmean <- .WeightedMean(data, lon = lon, lat = lat, mask = mask)
}
if(length(dim(data)) > 3) {
if (is.null(dim_names)) {
dim_names <- paste0("dim", 1:length(dim(data)))
......
man/ArrayToList.Rd 0 → 100644
View file @ dca0f085
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/ArrayToList.R
\name{ArrayToList}
\alias{ArrayToList}
\title{Split an array into list by a given array dimension}
\usage{
ArrayToList(data, dim, level = "list", names = NULL)
}
\arguments{
\item{data}{A multidimensional array.}
\item{dim}{A character string indicating the name of the dimension to split or an integer indicating the position of the dimension.}
\item{level}{A string character 'list' or 'sublist' indicating if it should be a list or a sublist. By default it creates a list.}
\item{names}{A vector of character strings to name the list (if it is a single string, it would be reused) or a single character string to name the elements in the sublist.}
}
\value{
A list of arrays of the length of the dimension set in parameter 'dim'.
}
\description{
This function splits an array into a list as required by PlotLayout function from package "s2dv" when parameter 'special_args' is used. The function ArrayToList allows to add names to the elements of the list in two different levels, the 'list' or the 'sublist'.
}
\examples{
data <- array(1:240, c(month = 12, member = 5, time = 4))
# Create a list:
datalist <- ArrayToList(data, dim = 'month', level = 'list', names = month.name)
class(datalist)
class(datalist[[1]])
str(datalist)
# Create a sublist:
datalist <- ArrayToList(data, dim = 'month', level = 'sublist', names = 'dots')
class(datalist)
class(datalist[[1]])
class(datalist[[1]][[1]])
str(datalist)
}
\seealso{
\link[s2dv]{PlotLayout}
}
vignettes/anomaly_agreement.Rmd
View file @ dca0f085
......@@ -32,7 +32,7 @@ All the other R packages involved can be installed directly from CRAN and loaded
```r
library(abind)
library(s2dverification)
library(s2dv)
library(ggplot2)
```
......@@ -170,18 +170,18 @@ model time lon lat
4 90 12 8
```
+ Now, the mean climatology value for each grid point and model can be computed by using the `Mean1Dim()` function belonging to the **s2dverification package**. The position of the temporal dimension should be specified in parameter `posdim`.
+ Now, the mean climatology value for each grid point and model can be computed by using the `MeanDims()` function belonging to the **s2dv package**. The position of the temporal dimension should be specified in parameter `posdim`.
```r
climatology <- Mean1Dim(summer_historical$data, posdim = 2)
climatology <- MeanDims(summer_historical$data, dims = 'time')
```
+ `Season()` function from **s2dverification package** returns the mean annual time series for the selected season by defining the parameters of the initial month of the data (`monini = 1`), the first month of the season (`moninf = 6`) and the final month of the season (`monsup = 8`). The last two parameters, `moninf` and `monsup`, have their origin with respect to the first one `monini`.
+ `Season()` function from **s2dv package** returns the mean annual time series for the selected season by defining the parameters of the initial month of the data (`monini = 1`), the first month of the season (`moninf = 6`) and the final month of the season (`monsup = 8`). The last two parameters, `moninf` and `monsup`, have their origin with respect to the first one `monini`.
```r
summer_projection <- Season(multimodel_projection, posdim = 3,
summer_projection <- Season(multimodel_projection, time_dim = 'time',
monini = 1, moninf = 6, monsup = 8)
```
......@@ -192,15 +192,17 @@ By running the next lines, it is possible to check the dimensions of the data:
model lon lat
4 12 8
> dim(summer_projection)
model var time lon lat
4 1 95 12 8
time model var lon lat
95 4 1 12 8
```
+ A new dimension will be added by running the function `InsertDim()` in order to obtain the same dimensions as the projections data. `InsertDim()` repeats the original data the required number of times (21 years of future simulations) in the adequated position (the temporal dimension in the summer_projection data is in the third position).
```r
climatology <- InsertDim(InsertDim(climatology, posdim = 2, lendim = 1),
posdim = 3, lendim = 95)
climatology <- InsertDim(
InsertDim(
climatology, posdim = 2, lendim = 1, name = 'var'),
posdim = 1, lendim = 95, name = 'time')
```
+ The anomaly for each model is obtained by simply subtracting.
......@@ -214,7 +216,7 @@ anomaly <- summer_projection - climatology
In order to obtain a spatial visualitzation, the temporal mean is computed. So, the time average anomalies for all models is saved in the `average` object. `AnoAgree()` function from **ClimProjDiags package** calculates the percentages of models which agrees with a positive or negative mean in each grid point.
```r
average <- Mean1Dim(anomaly, which(names(dim(anomaly)) == "time"))
average <- MeanDims(anomaly, dims = "time")
agreement <- AnoAgree(average, membersdim = which(names(dim(average)) == "model"))
```
......@@ -227,7 +229,7 @@ agreement_threshold <- 80
colorbar_lim <- ceiling(max(abs(max(average)), abs(min(average))))
brks <- seq(-colorbar_lim, colorbar_lim, length.out = 11)
PlotEquiMap(drop(Mean1Dim(average, which(names(dim(average)) == "model"))),
PlotEquiMap(drop(MeanDims(average, dims = "model")),
lat = lat, lon = lon, units = "K", brks = brks,
toptitle = paste(var, "- climatology:", start_climatology, "to",
end_climatology, "and future simulation:",
......@@ -243,17 +245,17 @@ PlotEquiMap(drop(Mean1Dim(average, which(names(dim(average)) == "model"))),
### 5- Multi-model agreement temporal visualization
To visualize the time evolution of multi-model agreement, the spatial average is performed by a grid pixel size using the `WeightedMean` function from the **ClimProjDiags package**. Also, a smooth filter is applied with the `Smoothing()` function from the **s2dverifiction package**. In this example, a 5-year moving window filter is applied by defining the parameter `runmeanlen = 5`.
To visualize the time evolution of multi-model agreement, the spatial average is performed by a grid pixel size using the `WeightedMean` function from the **ClimProjDiags package**. Also, a smooth filter is applied with the `Smoothing()` function from the **s2dv package**. In this example, a 5-year moving window filter is applied by defining the parameter `runmeanlen = 5`.
```r
temporal <- drop(WeightedMean(anomaly, lon = lon, lat = lat, mask = NULL))
temporal <- Smoothing(temporal, runmeanlen = 5, numdimt = 2)
temporal <- Smoothing(temporal, time_dim = 'time', runmeanlen = 5)
```
Before visualizing, a data frame with the proper format is created.
```r
data_frame <- as.data.frame.table(t(temporal))
data_frame <- as.data.frame.table(temporal)
years <- rep(start_projection : end_projection, 4)
data_frame$Year <- c(years)
names(data_frame)[2] <- "Model"
......
vignettes/diurnaltemp.Rmd
View file @ dca0f085
......@@ -32,7 +32,7 @@ library(ClimProjDiags)
All the other R packages involved can be installed directly from CRAN and loaded as follows:
```r
library(s2dverification)
library(s2dv)
library(abind)
library(parallel)
```
......@@ -78,12 +78,16 @@ for (j in 1 : 8) {
tmin_historical <- abind(tmin_historical, gridnew2, along = 3)
}
tmax_historical <- InsertDim(InsertDim(tmax_historical, posdim = 1, lendim = 1),
posdim = 1, lendim = 1)
tmin_historical <- InsertDim(InsertDim(tmin_historical, posdim = 1, lendim = 1),
posdim = 1, lendim = 1)
names(dim(tmax_historical)) <- c("model", "var", "time", "lon", "lat")
names(dim(tmin_historical)) <- c("model", "var", "time", "lon", "lat")
names(dim(tmax_historical)) <- c('time', 'lon', 'lat')
names(dim(tmin_historical)) <- c('time', 'lon', 'lat')
tmax_historical <- InsertDim(InsertDim(tmax_historical, posdim = 1,
lendim = 1, name = 'var'),
posdim = 1, lendim = 1, name = 'model')
tmin_historical <- InsertDim(InsertDim(tmin_historical, posdim = 1,
lendim = 1, name = 'var'),
posdim = 1, lendim = 1, name = 'model')
time <- seq(ISOdate(1971, 1, 1), ISOdate(2000, 12, 31), "day")
metadata <- list(time = list(standard_name = 'time', long_name = 'time',
calendar = 'proleptic_gregorian',
......@@ -125,12 +129,16 @@ for (j in 1 : 8) {
tmin_projection <- abind(tmin_projection, gridnew2, along = 3)
}
tmax_projection <- InsertDim(InsertDim(tmax_projection, posdim = 1, lendim = 1),
posdim = 1, lendim = 1)
tmin_projection <- InsertDim(InsertDim(tmin_projection, posdim = 1, lendim = 1),
posdim = 1, lendim = 1)
names(dim(tmax_projection)) <- c("model", "var", "time", "lon", "lat")
names(dim(tmin_projection)) <- c("model", "var", "time", "lon", "lat")
names(dim(tmax_projection)) <- c("time", "lon", "lat")
names(dim(tmin_projection)) <- c("time", "lon", "lat")
tmax_projection <- InsertDim(InsertDim(tmax_projection, posdim = 1,
lendim = 1, name = 'var'),
posdim = 1, lendim = 1, name = 'model')
tmin_projection <- InsertDim(InsertDim(tmin_projection, posdim = 1,
lendim = 1, name = 'var'),
posdim = 1, lendim = 1, name = 'model')
time <- seq(ISOdate(2006, 1, 1), ISOdate(2100, 12, 31), "day")
metadata <- list(time = list(standard_name = 'time', long_name = 'time',
calendar = 'proleptic_gregorian',
......@@ -193,7 +201,7 @@ A four panel plot can be generated to visualize the seasonal indicator by runnin
```r
dtr_rcp <- array(dim = c(length(lon), length(lat), 4))
for (j in 1:4){
dtr_rcp[,,j] <- Mean1Dim(dtr_indicator$indicator[,j,,,,], 1)
dtr_rcp[,,j] <- MeanDims(dtr_indicator$indicator[,j,,,,], 'year')
}
breaks <- 0 : (max(dtr_rcp) + 1)
......@@ -224,10 +232,10 @@ dtr_indicator_reference <- DTRIndicator(tmax = tmax_historical,
The comparison between the reference and the future projection diurnal temperature variation indicator can be computed With a simple subtraction. To visualize the result, `PlotLayout` function will be applied again. The resulting plot will be saved in the working directory.
```r
dtr_diff <- array(dim=c(length(lat), length(lon), 4))
dtr_diff <- array(dim = c(length(lat), length(lon), 4))
for (i in 1:4){
dtr_diff[,,i] <- Mean1Dim(dtr_indicator$indicator[,i,1,1,,], 1) -
Mean1Dim(dtr_indicator_reference$indicator[,i,1,1,,], 1)
dtr_diff[,,i] <- MeanDims(dtr_indicator$indicator[, i, 1, 1, , ], 'year') -
MeanDims(dtr_indicator_reference$indicator[, i, 1, 1, , ], 'year')
}
PlotLayout(PlotEquiMap, c(1, 2), lon = lon, lat = lat, var = dtr_diff,
......
vignettes/extreme_indices.Rmd
View file @ dca0f085
......@@ -39,7 +39,7 @@ All the other R packages involved can be installed directly from CRAN and loaded
```r
library(s2dverification)
library(s2dv)
library(abind)
library(multiApply)
library(ggplot2)
......@@ -71,9 +71,11 @@ for (j in 1 : 8) {
gridnew <- apply(gridlon, 2, function(x) {x - rnorm(10958, mean = j * 0.5, sd = 3)})
tmax_historical <- abind(tmax_historical, gridnew, along = 3)
}
tmax_historical <- InsertDim(InsertDim(tmax_historical, posdim = 1, lendim = 1),
posdim = 1, lendim = 1)
names(dim(tmax_historical)) <- c("model", "var", "time", "lon", "lat")
names(dim(tmax_historical)) <- c("time", "lon", "lat")
tmax_historical <- InsertDim(InsertDim(tmax_historical, posdim = 1,
lendim = 1, name = 'var'),
posdim = 1, lendim = 1, name = 'model')
time <- seq(ISOdate(1971, 1, 1), ISOdate(2000, 12, 31), "day")
metadata <- list(time = list(standard_name = 'time', long_name = 'time',
calendar = 'proleptic_gregorian',
......@@ -102,8 +104,10 @@ for (j in 1 : 8) {
rnorm(34698, mean = j * 0.5, sd = 3)})
tmax_projection <- abind(tmax_projection, gridnew, along = 3)
}
tmax_projection <- InsertDim(InsertDim(tmax_projection, posdim = 1, lendim = 1), posdim = 1, lendim = 1)
names(dim(tmax_projection)) <- c("model", "var", "time", "lon", "lat")
names(dim(tmax_projection)) <- c("time", "lon", "lat")
tmax_projection <- InsertDim(InsertDim(tmax_projection, posdim = 1,
lendim = 1, name = 'var'),
posdim = 1, lendim = 1, name = 'model')
time <- seq(ISOdate(2006, 1, 1), ISOdate(2100, 12, 31), "day")
metadata <- list(time = list(standard_name = 'time', long_name = 'time',
calendar = 'proleptic_gregorian',
......@@ -159,12 +163,12 @@ names(dim(anomaly_data))[1] <- "time"
```
This data can be detrended by applying the `Trend` function from **s2dverification package**. In order to remove the trend from the `tmax_historical`, the correction is calculated by subtracting the `detrended_data` to the `anomaly_data`.
This data can be detrended by applying the `Trend` function from **s2dv package**. In order to remove the trend from the `tmax_historical`, the correction is calculated by subtracting the `detrended_data` to the `anomaly_data`.
```r
detrended_data <- Trend(anomaly_data,
posTR = which(names(dim(anomaly_data)) == "time"))
detrended_data <- Trend(anomaly_data, time_dim = "time")
diff <- anomaly_data - detrended_data$detrended
diff <- aperm(diff, c(2,3,1,4,5))
detrended_data <- tmax_historical - diff
......@@ -247,8 +251,7 @@ masc <- rep(0, 8 * 12)
masc[c(5 : 12, 18 : 24, 31 : 34, 43, 44, 47, 56 : 60, 67 : 72, 79,
82 : 84, 93 : 96)] <- 1
dim(masc) <- c(12, 8)
PlotEquiMap(Mean1Dim(HeatExtremeIndex,
which(names(dim(HeatExtremeIndex)) == "year")),
PlotEquiMap(MeanDims(HeatExtremeIndex, "year"),
lon = lon, lat = lat, filled.continents = FALSE,
toptitle = "Extreme Heat Index", dots = masc,
fileout = "SpatialExtremeHeatIndex.png")
......@@ -263,7 +266,7 @@ The inland average of the Extreme Heat Index can be computed to plot its time ev
```r
temporal <- WeightedMean(HeatExtremeIndex, lon = lon, lat = lat, mask = drop(masc))
temporal_3ysmooth <- Smoothing(temporal, runmeanlen = 3, numdimt = 1)
temporal_3ysmooth <- Smoothing(temporal, runmeanlen = 3, time_dim = 'year')
```
......@@ -335,8 +338,7 @@ Spatial representation of the Extreme Drought Index:
```r
PlotEquiMap(Mean1Dim(DroughtExtremeIndex,
which(names(dim(DroughtExtremeIndex)) == "year")),
PlotEquiMap(MeanDims(DroughtExtremeIndex, "year"),
lon = lon, lat = lat, filled.continents = FALSE,
toptitle = "Drought Index", brks = seq(-1, 1, 0.01),
fileout = "SpatialDroughtIndex.png")
......@@ -352,7 +354,7 @@ Evolution of inland average of the Extreme Drought Index:
```r
temporal <- WeightedMean(DroughtExtremeIndex, lon = lon, lat = lat,
mask = drop(masc))
temporal_5ysmooth <- Smoothing(temporal, runmeanlen = 5, numdimt = 1)
temporal_5ysmooth <- Smoothing(temporal, runmeanlen = 5, time_dim = "year")
png("Temporal_Inland_ExtremeDroughtIndex.png", width = 8, height = 5, units= 'in',
res = 100, type = "cairo")
plot(2006: 2100, temporal, type = "l", lty = 5, lwd = 2, bty = 'n',
......@@ -402,9 +404,8 @@ Spatial representation of the Extreme Flooding Index:
```r
PlotEquiMap(Mean1Dim(FloodingExtremeIndex,
which(names(dim(FloodingExtremeIndex)) == "year")), lon = lon,
lat = lat, filled.continents = FALSE,
PlotEquiMap(MeanDims(FloodingExtremeIndex, "year"),
lon = lon, lat = lat, filled.continents = FALSE,
toptitle = "Extreme Flooding Index",
brks = seq(-1, 1, 0.1), fileout = "SpatialFloodingIndex.png")
```
......@@ -419,7 +420,7 @@ PlotEquiMap(Mean1Dim(FloodingExtremeIndex,
```r
temporal <- WeightedMean(FloodingExtremeIndex, lon = lon, lat = lat,
mask = drop(masc))
temporal_3ysmooth <- Smoothing(temporal, runmeanlen = 3, numdimt = 1)
temporal_3ysmooth <- Smoothing(temporal, runmeanlen = 3, time_dim = "year")
png("Temporal_Inland_ExtremeFloodingIndex.png", width = 8, height = 5,
units= 'in', res = 100, type = "cairo")
plot(2006 : 2100, temporal, type = "l", lty = 5, lwd = 2, bty = 'n',
......@@ -460,7 +461,7 @@ A spatial visualization can be performed by executing:
```r
PlotEquiMap(Mean1Dim(aci, which(names(dim(aci)) == "year")), lon = lon,
PlotEquiMap(MeanDims(aci, "year"), lon = lon,
lat = lat, filled.continents = FALSE, toptitle = "Indices Combination",
fileout = "CombinedIndices.png")
```
......
vignettes/heatcoldwaves.Rmd
View file @ dca0f085
......@@ -32,7 +32,7 @@ All the other R packages involved can be installed directly from CRAN and loaded
```r
library(abind)
library(s2dverification)
library(s2dv)
library(parallel)
```
......@@ -68,9 +68,10 @@ for (j in 1 : 8) {
gridnew <- apply(gridlon, 2, function(x) {x - rnorm(10958, mean = j * 0.5, sd = 3)})
tmax_historical <- abind(tmax_historical, gridnew, along = 3)
}
tmax_historical <- InsertDim(InsertDim(tmax_historical, posdim = 1, lendim = 1),
posdim = 1, lendim = 1)
names(dim(tmax_historical)) <- c("model", "var", "time", "lon", "lat")
names(dim(tmax_historical)) <- c("time", "lon", "lat")
tmax_historical <- InsertDim(InsertDim(tmax_historical, posdim = 1,
lendim = 1, name = "var"),
posdim = 1, lendim = 1, name = "model")
time <- seq(ISOdate(1971, 1, 1), ISOdate(2000, 12, 31), "day")
metadata <- list(time = list(standard_name = 'time', long_name = 'time',
calendar = 'proleptic_gregorian',
......@@ -97,9 +98,10 @@ for (j in 1 : 8) {
sd = 3)})
tmax_projection <- abind(tmax_projection, gridnew, along = 3)
}
tmax_projection <- InsertDim(InsertDim(tmax_projection, posdim = 1, lendim = 1),
posdim = 1, lendim = 1)
names(dim(tmax_projection)) <- c("model", "var", "time", "lon", "lat")
names(dim(tmax_projection)) <- c("time", "lon", "lat")
tmax_projection <- InsertDim(InsertDim(tmax_projection, posdim = 1,
lendim = 1, name = "var"),
posdim = 1, lendim = 1, name = "model")
time <- seq(ISOdate(2006, 1, 1), ISOdate(2100, 12, 31), "day")
metadata <- list(time = list(standard_name = 'time', long_name = 'time',
calendar = 'proleptic_gregorian',
......