DESCRIPTION

 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

 # Generated by roxygen2: do not edit by hand
 
 export(AnoAgree)
+export(ArrayToList)
 export(Climdex)
 export(CombineIndices)
 export(DTRIndicator)

R/ArrayToList.R

+#' 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

@@ -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

+% 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

@@ -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

@@ -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

@@ -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

@@ -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',