Merge branch 'devel' of github.com:SantanderMetGroup/transformeR into…

… devel

DESCRIPTION

@@ -23,8 +23,8 @@ Suggests:
     visualizeR
 Type: Package
 Title: A climate4R package for general climate data manipulation and transformation
-Version: 2.0.2
-Date: 2020-09-02
+Version: 2.1.2
+Date: 2021-07-07
 Authors@R: c(person("Bedia", "Joaquín", email = "[email protected]", role = c("aut","cre"), comment = c(ORCID = "0000-0001-6219-4312")),
              person("Jorge", "Baño Medina", email = "[email protected]", role = "aut"),
              person("Ana", "Casanueva", email = "[email protected]", role = "aut"),

NAMESPACE

@@ -21,7 +21,9 @@ export(dataSplit)
 export(detrendGrid)
 export(easyVeri.signif)
 export(easyVeri2grid)
+export(fillGrid)
 export(fillGridDates)
+export(fillGridSpatial)
 export(filterGrid)
 export(filterNA)
 export(get2DmatCoordinates)

NEWS

@@ -120,9 +120,20 @@ This release of transformeR does not ensure backwards compatibility. Some deprec
  * Force nearest-neighbour application among non-overlapping domains in `interpGrid`
  * Other minor documentation updates
 
-## 2.0.2 (in devel)
+## 2.1.0 (17 Mar 2021)
+ * Implementation of Lamb Weather Typing for southern hemisphere
+ * Inclusion of type 27 (unclassified) in Lamb WT
+ * New functionality for spatial aggregation of rectangular domains
+ * Bug fix in partially masked latitudinal weighted domain averages.
+ * New option to use masking in interpGrid
  * Bug fix in subsetVar attributes
- * Internal falg in bindGrid.member to skip temporal checks
- * Minor changes and documentation updates
+ * Internal flag in bindGrid.member to skip temporal checks
+ * Other minor changes and documentation updates
 
+## 2.1.1 (31 May 2021)
+ * New dateutils-family helper functions for consistent date-time class handling (including time-zone defs)
+ * Consistent treatment of multiple time resolutions in fillGridDates 
 
+## 2.1.2 (07 Jul 2021)
+ * Allow for non-rescaled EOF recovery in EOF2clim
+ * Ensure regular grids in upscaleGrid outputs 

R/EOF2clim.R

@@ -2,14 +2,29 @@
 #' @description Converts a \code{prinComp} output to a climatology of EOFs. 
 #' @param prinCompObj PCA object
 #' @param ind.var index position of the target variable
-#' @param member index position of the member whose EOFs are to be displayed
-#' @param n.eofs Number of EOFs to be displayed (from 1 to n.eofs)
+#' @param member Integer. index position of the member whose EOFs are to be displayed.
+#'  Default to 1 (ignored if no members available)
+#' @param n.eofs Integer. Number of EOFs to be displayed (from 1 to \code{n.eofs}). Default to one -first- EOF.
+#' @param rescale Logical flag. Default to \code{TRUE}, which rescales the EOF to the original input variable units.
 #' @return A climatological grid. Note that EOFs are treated as members
 #' @importFrom magrittr %>% 
 #' @export
 #' @author J Bedia
+#' @examples \dontrun{
+#' require(climate4R.datasets)
+#' data("NCEP_Iberia_hus850", "NCEP_Iberia_psl", "NCEP_Iberia_ta850")
+#' multigrid <- makeMultiGrid(NCEP_Iberia_hus850, NCEP_Iberia_psl, NCEP_Iberia_ta850)
+#' # In this example, we retain the PCs explaining the 99\% of the variance
+#' pca <- prinComp(multigrid, v.exp = c(.95,0.90,.90), keep.orig = FALSE)
+#' require(visualizeR)
+#' # Original EOF
+#' EOF2clim(pca, ind.var = "psl") %>% spatialPlot(backdrop.theme = "coastline")
+#' # Rescaled EOF
+#' EOF2clim(pca, ind.var = "psl", rescale = FALSE) %>% spatialPlot(backdrop.theme = "coastline")
+#' }
 
-EOF2clim <- function(prinCompObj, ind.var, member, n.eofs) {
+
+EOF2clim <- function(prinCompObj, ind.var, member = 1L, n.eofs = 1L, rescale = TRUE) {
   varNames <- attributes(prinCompObj)$names
   levs <- attributes(prinCompObj)$level
   x <- attributes(prinCompObj)$xCoords
@@ -22,7 +37,10 @@ EOF2clim <- function(prinCompObj, ind.var, member, n.eofs) {
   eofs <- prinCompObj[[ind.var]][[member]]$EOFs[, 1:n.eofs, drop = FALSE]
   prinCompObj <- NULL
   # Rescale EOFs
-  aux <- (eofs * sigma + mu) %>% t() %>% mat2Dto3Darray(x, y) %>% list()
+  if (isTRUE(rescale)) {
+    eofs <- (eofs * sigma + mu)
+  }
+  aux <- t(eofs) %>% mat2Dto3Darray(x, y) %>% list()
   # Recover grids structure (EOFS are treated as members, time = 1 like a climatology
   Data <- do.call("abind", c(aux, along = -1)) %>% unname()
   attr(Data, "dimensions") <- c("time", "member", "lat", "lon")

R/aggregateGrid.R

@@ -1,6 +1,6 @@
 #     aggregateGrid.R Flexible grid aggregation along selected dimensions
 #
-#     Copyright (C) 2017 Santander Meteorology Group (http://www.meteo.unican.es)
+#     Copyright (C) 2021 Santander Meteorology Group (http://www.meteo.unican.es)
 #
 #     This program is free software: you can redistribute it and/or modify
 #     it under the terms of the GNU General Public License as published by
@@ -20,16 +20,19 @@
 #' @description Aggregates a grid along the target dimensions using user-defined functions.
 #' @param grid a grid or multigrid to be aggregated.
 #' @param aggr.d Daily aggregation function (for sub-daily data only). A list indicating the name of the
-#'  aggregation function in first place, and other optional arguments to be passed to the aggregation function. See the examples.
+#'  aggregation function in first place, and other optional arguments to be passed to the aggregation function.
+#'  To be on the safe side, the function in \code{FUN} should be always indicated as a character string. See the examples.
 #' @param aggr.m Same as \code{aggr.d}, but indicating the monthly aggregation function. 
 #' @param aggr.y Same as \code{aggr.d}, but indicating the annual aggregation function. 
 #' @param aggr.s Same as \code{aggr.d}, but indicating the seasonal aggregation function. The season can be indicated
 #' as shown in this example:  aggr.s = list(FUN = list("mean", na.rm = TRUE), season = c(12,1,2)) 
 #' @param aggr.mem Same as \code{aggr.d}, but indicating the function for computing the member aggregation.
-#' @param aggr.lat Same as \code{aggr.d}, indicating the aggregation function to be applied along latitude.
+#' @param aggr.spatial Same as \code{aggr.d}, but indicating the aggregation function in case of rectangular domains to be aggregated
+#' as a unique time series grid (or multimember time series grid.)
+#' @param aggr.lat Same as \code{aggr.d}, indicating the aggregation function to be applied along latitude only.
 #' @param weight.by.lat Logical. Should latitudinal averages be weighted by the cosine of latitude?.
-#' Default to \code{TRUE}. Ignored if no \code{aggr.lat} function is indicated, or a function different from \code{"mean"}
-#' is applied.
+#' Default to \code{TRUE}. Ignored if no \code{aggr.lat} or \code{aggr.spatial} function is indicated,
+#'  or a function different from \code{"mean"} is applied.
 #' @param aggr.lon Same as \code{aggr.lat}, but for longitude.
 #' @param aggr.loc Same as \code{aggr.d}, indicating the aggregation function to be applied along the loc dimension.
 #' @template templateParallelParams
@@ -95,6 +98,7 @@ aggregateGrid <- function(grid,
                           aggr.m = list(FUN = NULL),
                           aggr.y = list(FUN = NULL),
                           aggr.s = list(FUN = NULL, season = NULL),
+                          aggr.spatial = list(FUN = NULL),
                           aggr.lat = list(FUN = NULL),
                           weight.by.lat = TRUE,
                           aggr.lon = list(FUN = NULL),
@@ -124,6 +128,11 @@ aggregateGrid <- function(grid,
                 climatology(clim.fun = aggr.s$FUN)  
         }) %>% bindGrid(dimension = "time")
     }
+    if (!is.null(aggr.spatial$FUN)) {
+        grid <- spatialAggregation(grid, aggr.fun = aggr.spatial,
+                                   weight.by.lat = weight.by.lat,
+                                   parallel, max.ncores, ncores)
+    }
     if (!is.null(aggr.lat$FUN)) {
         grid <- latAggregation(grid, aggr.lat, weight.by.lat, parallel, max.ncores, ncores)
     }
@@ -275,6 +284,7 @@ timeAggregation <- function(grid, aggr.type = c("DD","MM","YY"), aggr.fun, paral
 
 
 latAggregation <- function(grid, aggr.fun, weight.by.lat, parallel, max.ncores, ncores) {
+    stopifnot(is.logical(weight.by.lat))
     dimNames <- getDim(grid)
     if (!"lat" %in% dimNames) {
         message("There is not lat dimension: 'aggr.lat' option was ignored.")
@@ -380,3 +390,48 @@ locAggregation <- function(grid, aggr.fun, weight.by.lat, parallel, max.ncores,
       }
       return(grid)
 }
+
+#' @title Spatial Aggregation
+#' @description Spatial aggregation for rectangular domains. 
+#' @param grid Input grid. 
+#' @param aggr.fun Aggregation function
+#' @param weight.by.lat Logical flag
+#' @return A spatially averaged time series grid (possibly multimember)
+#' @keywords internal
+#' @author J Bedia
+
+spatialAggregation <- function(grid, aggr.fun, weight.by.lat, parallel, max.ncores, ncores) {
+  stopifnot(is.logical(weight.by.lat))
+  dimNames <- getDim(grid)
+  if (!any(c("lon", "lat") %in% dimNames)) {
+    message("\'lon\' and/or \'lat\' dimensions are missing in the input grid: 'aggr.spatial' option was ignored.")
+  } else {
+    if (isTRUE(weight.by.lat)) {
+      lat.weights <- latWeighting(grid)
+      weight.matrix <- matrix(lat.weights, nrow = getShape(grid, "lat"), ncol = getShape(grid, "lon"))
+      if (aggr.fun[["FUN"]] == "mean") {
+        message("Calculating areal weights...")
+        aggr.fun <- list(FUN = "weighted.mean", w = weight.matrix, na.rm = TRUE)
+      } else {
+        message("Spatial weighting skipped: It only applies to \'mean\' aggregation function")
+      }
+    }
+    parallel.pars <- parallelCheck(parallel, max.ncores, ncores)
+    mar <- grep("lat|lon", dimNames, invert = TRUE)
+    aggr.fun[["MARGIN"]] <- mar
+    aggr.fun[["X"]] <- grid$Data
+    out <- if (parallel.pars$hasparallel) {
+      message("[", Sys.time(), "] - Aggregating lat dimension in parallel...")
+      on.exit(parallel::stopCluster(parallel.pars$cl))
+      aggr.fun[["cl"]] <- parallel.pars$cl
+      do.call("parApply", aggr.fun)
+    } else {
+      message("[", Sys.time(), "] - Aggregating spatially...")
+      do.call("apply", aggr.fun)
+    }
+    grid$Data <- out
+    attr(grid$Data, "dimensions") <- dimNames[mar]
+    message("[", Sys.time(), "] - Done.")
+  }
+  return(grid)
+}

R/binaryGrid.R

@@ -72,7 +72,7 @@ binaryGrid <- function(x,
         if (isRegular(x)) {
           xx <- suppressWarnings(array3Dto2Dmat(redim(subsetGrid(x,members = j), member = FALSE)$Data))
         } else {
-          xx <- x$Data[j,,]
+          xx <- x$Data[j,,] %>% as.matrix()
         }
         s <- matrix(runif(nrow(xx)*ncol(xx),min = 0,max = 1),nrow = nrow(xx), ncol = ncol(xx))
         xbin <- (xx > s)*1

R/bindGrid.R

@@ -1,6 +1,6 @@
 #     bindGrid.R Grid binding along user-defined dimension
 #
-#     Copyright (C) 2020 Santander Meteorology Group (http://www.meteo.unican.es)
+#     Copyright (C) 2021 Santander Meteorology Group (http://www.meteo.unican.es)
 #
 #     This program is free software: you can redistribute it and/or modify
 #     it under the terms of the GNU General Public License as published by
@@ -354,51 +354,51 @@ sortDim.spatial <- function(grid, dimension = c("y", "x")) {
 #' @author M De Felice, J Bedia
 
 bindGrid.time <- function(..., tol) {
-  grid.list <- list(...)
-  if (length(grid.list) == 1) {
-    grid.list <- unlist(grid.list, recursive = FALSE)
-  }
-  if (length(grid.list) < 2) {
-    stop("The input must be a list of at least two grids")
-  }
-  grid.list <- lapply(grid.list, function(i){
-    loc <- !isRegular(i)
-    redim(i, loc = loc, var = TRUE)
-  })
-    for (i in 2:length(grid.list)) {
+   grid.list <- list(...)
+   if (length(grid.list) == 1) {
+      grid.list <- unlist(grid.list, recursive = FALSE)
+   }
+   if (length(grid.list) < 2) {
+      stop("The input must be a list of at least two grids")
+   }
+   grid.list <- lapply(grid.list, function(i){
+      loc <- !isRegular(i)
+      redim(i, loc = loc, var = TRUE)
+   })
+   for (i in 2:length(grid.list)) {
       # Spatial test
       if (!isTRUE(all.equal(grid.list[[1]]$xyCoords, grid.list[[i]]$xyCoords, check.attributes = FALSE, tolerance = tol))) {
-        stop("Input data is not spatially consistent")
+         stop("Input data is not spatially consistent")
       }
       # Member
       if (getShape(grid.list[[1]])[match('member', getDim(grid.list[[1]]))] != getShape(grid.list[[i]])[match('member', getDim(grid.list[[i]]))]) {
-        stop("Member dimension is not spatially consistent")
+         stop("Member dimension is not spatially consistent")
       }
-    }
-    ref <- grid.list[[1]]
-    dimNames <- getDim(ref) 
-    dim.bind <- grep("^time", dimNames)
-    data.list <- lapply(grid.list, FUN = "[[", "Data")
-    ref[["Data"]] <- unname(do.call("abind", c(data.list, along = dim.bind)))
-    data.list <- NULL
-    start.list <- lapply(grid.list, FUN = function(x) {
+   }
+   ref <- grid.list[[1]]
+   dimNames <- getDim(ref) 
+   dim.bind <- grep("^time", dimNames)
+   data.list <- lapply(grid.list, FUN = "[[", "Data")
+   ref[["Data"]] <- unname(do.call("abind", c(data.list, along = dim.bind)))
+   data.list <- NULL
+   start.list <- lapply(grid.list, FUN = function(x) {
       getRefDates(x)
-    })
-    end.list <- lapply(grid.list, FUN = function(x) {
+   })
+   end.list <- lapply(grid.list, FUN = function(x) {
       getRefDates(x, "end")
-    })
-    grid.list <- NULL
-    refdates <- list(start = do.call(c, start.list),
-                     end = do.call(c, end.list))
-    attr(ref[["Data"]], "dimensions") <- dimNames
-    n.vars <- getShape(ref, "var")
-    if (n.vars > 1) refdates <- rep(list(refdates), n.vars)
-    ref[["Dates"]] <- refdates
-    ref <- tryCatch({sortDim.time(ref)}, error = function(err) {
+   })
+   grid.list <- NULL
+   refdates <- list(start = do.call("c", start.list),
+                    end = do.call("c", end.list))
+   attr(ref[["Data"]], "dimensions") <- dimNames
+   n.vars <- getShape(ref, "var")
+   if (n.vars > 1) refdates <- rep(list(refdates), n.vars)
+   ref[["Dates"]] <- refdates
+   ref <- tryCatch({sortDim.time(ref)}, error = function(err) {
       warning("time dimension could not be sorted!")
       ref})
-    ref <- redim(ref, drop = TRUE)
-    return(ref)
+   ref <- redim(ref, drop = TRUE)
+   return(ref)
 }
 #end