diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R index 77189e4a2436426ca2eed3d8314610babd4a0fae..a7dd5ff3b592498452cd17f74d08234fb0da35d2 100644 --- a/R/CST_RainFARM.R +++ b/R/CST_RainFARM.R @@ -37,7 +37,8 @@ #' If omitted one of c("ftime", "sdate", "time") is searched and the first one with more #' than one element is chosen. #' @param verbose Logical for verbose output (default: FALSE). -#' @param drop_realization_dim Logical to remove the "realization" stochastic ensemble dimension (default: FALSE) +#' @param drop_realization_dim Logical to remove the "realization" stochastic ensemble dimension, +#' needed for saving data through function CST_SaveData (default: FALSE) #' with the following behaviour if set to TRUE: #' #' 1) if \code{nens==1}: the dimension is dropped; @@ -46,7 +47,8 @@ #' the "realization" and "member" dimensions are compacted (multiplied) and the resulting dimension is named "member"; #' #' 3) if \code{nens>1} and a "member" dimension does not exist: the "realization" dimension is renamed to "member". -#' +#' @param nprocs The number of parallel processes to spawn for the use for parallel computation in multiple cores. (default: 1) +#' #' @return CST_RainFARM() returns a downscaled CSTools object (i.e., of the #' class 's2dv_cube'). #' If \code{nens>1} an additional dimension named "realization" is added to the @@ -56,7 +58,7 @@ #' @import multiApply #' @import rainfarmr #' @examples -#' #Example using CST_RainFARM for a CSTools object +#' #Example 1: using CST_RainFARM for a CSTools object #' nf <- 8 # Choose a downscaling by factor 8 #' exp <- 1 : (2 * 3 * 4 * 8 * 8) #' dim(exp) <- c(dataset = 1, member = 2, sdate = 3, ftime = 4, lat = 8, lon = 8) @@ -70,21 +72,35 @@ #' res <- CST_RainFARM(data, nf, ww, nens=3) #' str(res) #' #List of 3 -#' # $ data: num [1:4, 1:64, 1:64, 1, 1:2, 1:3] 201 115 119 307 146 ... +#' # $ data: num [1, 1:2, 1:3, 1:3, 1:4, 1:64, 1:64] 260 553 281 278 143 ... #' # $ lon : num [1:64] 9.78 9.84 9.91 9.97 10.03 ... #' # $ lat : num [1:64] 39.8 39.8 39.9 40 40 ... #' dim(res$data) -#' # dataset member sdate ftime lat lon realization -#' # 1 2 3 4 64 64 3 -#'@export +#' # dataset member realization sdate ftime lat lon +#' # 1 2 3 3 4 64 64 +#' +#' #Example 2: using CST_RainFARM for a CSTools object with parallel processing, +#' #dropping the "realization" dimension to be able to save results using +#' #\code{CST_SaveExp}. +#' #Load dataset included in CSTools pacakge +#' data <- lonlat_prec +#' nf <- 8 +#' # Create a test array of weights +#' ww <- array(1., dim = c(dim(data$lon) * nf, dim(data$lat) * nf)) +#' res <- CST_RainFARM(data, nf, weights=ww, nprocs=2, drop_realization_dim = TRUE) +#' dim(res$data) +#' # dataset member sdate ftime lat lon +#' # 1 6 3 31 32 32 + +#' @export CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1, nens = 1, fglob = FALSE, fsmooth = TRUE, - time_dim = NULL, verbose = FALSE, + nprocs = 1, time_dim = NULL, verbose = FALSE, drop_realization_dim = FALSE) { res <- RainFARM(data$data, data$lon, data$lat, nf, weights, nens, slope, kmin, fglob, fsmooth, - time_dim, lon_dim = "lon", lat_dim = "lat", + nprocs, time_dim, lon_dim = "lon", lat_dim = "lat", drop_realization_dim, verbose) data$data <- res$data @@ -135,7 +151,7 @@ CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1, #' @param lat_dim Name of lat dimension ("lat" by default). #' @param verbose logical for verbose output (default: FALSE). #' @param drop_realization_dim Logical to remove the "realization" stochastic ensemble dimension (default: FALSE) -# with the following behaviour if set to TRUE: +#' with the following behaviour if set to TRUE: #' #' 1) if \code{nens==1}: the dimension is dropped; #' @@ -144,6 +160,7 @@ CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1, #' #' 3) if \code{nens>1} and a "member" dimension does not exist: the "realization" dimension is renamed to "member". #' +#' @param nprocs The number of parallel processes to spawn for the use for parallel computation in multiple cores. (default: 1) #' @return RainFARM() returns a list containing the fine-scale longitudes, latitudes #' and the sequence of \code{nens} downscaled fields. #' If \code{nens>1} an additional dimension named "realization" is added to the output array @@ -184,7 +201,7 @@ CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1, #' RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1, slope = 0, kmin = 1, fglob = FALSE, fsmooth = TRUE, - time_dim = NULL, lon_dim = "lon", lat_dim = "lat", + nprocs = 1, time_dim = NULL, lon_dim = "lon", lat_dim = "lat", drop_realization_dim = FALSE, verbose = FALSE) { # Ensure input grid is square and with even dimensions @@ -251,7 +268,8 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1, # Repeatedly apply .RainFARM result <- Apply(data, c(lon_dim, lat_dim, "rainfarm_samples"), .RainFARM, weights, nf, nens, slope, kmin, - fglob, fsmooth, verbose)$output1 + fglob, fsmooth, ncores = nprocs, verbose, + split_factor = "greatest")$output1 # result has dims: lon, lat, rainfarm_samples, realization, other dims # Expand back rainfarm_samples to compacted dims dim(result) <- c(dim(result)[1:2], cdim[-ind], dim(result)[-(1:3)]) @@ -259,11 +277,14 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1, # Reorder as it was in original data # + realization dim after member if it exists ienspos <- which(names(cdim0) == "member") - if ( length(ienspos) == 0) ienspos <- length(names(cdim0)) + if (length(ienspos) == 0) ienspos <- length(names(cdim0)) iorder <- sapply(c(names(cdim0)[1:ienspos], "realization", names(cdim0)[-(1:ienspos)]), grep, names(dim(result))) + ndim <- names(dim(result)) result <- aperm(result, iorder) + # R < 3.2.3 compatibility fix + names(dim(result)) <- ndim[iorder] if (drop_realization_dim) { cdim <- dim(result) @@ -295,7 +316,7 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1, #' @param kmin First wavenumber for spectral slope (default: \code{kmin=1}). #' @param nens Number of ensemble members to produce (default: \code{nens=1}). #' @param fglob Logical to conseve global precipitation over the domain (default: FALSE). -#' @param fsmooth Logical to conserve precipitation with a smnoothing kernel (default: TRUE). +#' @param fsmooth Logical to conserve precipitation with a smoothing kernel (default: TRUE). #' @param verbose Logical for verbose output (default: FALSE). #' @return .RainFARM returns a downscaled array with dimensions (lon, lat, time, realization) #' @noRd @@ -324,15 +345,17 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1, # and array indexing. Derived from Stack Overflow issue # https://stackoverflow.com/questions/14500707/select-along-one-of-n-dimensions-in-array .subset <- function(field, dim_name, range, drop = FALSE) { - - idim <- which( names(dim(field)) %in% dim_name ) + + ndim <- names(dim(field)) + idim <- which(ndim %in% dim_name ) # Create list representing arguments supplied to [ # bquote() creates an object corresponding to a missing argument indices <- rep(list(bquote()), length(dim(field))) indices[[idim]] <- range - # do.call on the indices - field <- do.call("[",c(list(field), indices, list(drop = drop))) - + field <- do.call("[", c(list(field), indices, list(drop = drop))) + # Needed for R <=3.2 + names(dim(field)) <- ndim + return(field) } diff --git a/man/CST_RainFARM.Rd b/man/CST_RainFARM.Rd index ca32e2d8f0f37b36a434db4554dfdd1eb2e2675e..53af66ee23d211c1721d7c05a1091a510bd7aca5 100644 --- a/man/CST_RainFARM.Rd +++ b/man/CST_RainFARM.Rd @@ -5,8 +5,8 @@ \title{RainFARM stochastic precipitation downscaling of a CSTools object} \usage{ CST_RainFARM(data, nf, weights = 1, slope = 0, kmin = 1, nens = 1, - fglob = FALSE, fsmooth = TRUE, time_dim = NULL, verbose = FALSE, - drop_realization_dim = FALSE) + fglob = FALSE, fsmooth = TRUE, nprocs = 1, time_dim = NULL, + verbose = FALSE, drop_realization_dim = FALSE) } \arguments{ \item{data}{An object of the class 's2dv_cube' as returned by `CST_Load`, @@ -36,6 +36,8 @@ meaning that the slope is determined automatically over the dimensions specified \item{fsmooth}{Logical to conserve precipitation with a smoothing kernel (default: TRUE).} +\item{nprocs}{The number of parallel processes to spawn for the use for parallel computation in multiple cores. (default: 1)} + \item{time_dim}{String or character array with name(s) of dimension(s) (e.g. "ftime", "sdate", "member" ...) over which to compute spectral slopes. If a character array of dimension names is provided, the spectral slopes @@ -45,7 +47,8 @@ than one element is chosen.} \item{verbose}{Logical for verbose output (default: FALSE).} -\item{drop_realization_dim}{Logical to remove the "realization" stochastic ensemble dimension (default: FALSE) +\item{drop_realization_dim}{Logical to remove the "realization" stochastic ensemble dimension, +needed for saving data through function CST_SaveData (default: FALSE) with the following behaviour if set to TRUE: 1) if \code{nens==1}: the dimension is dropped; @@ -71,7 +74,7 @@ Adapted for climate downscaling and including orographic correction as described in Terzago et al. 2018. } \examples{ -#Example using CST_RainFARM for a CSTools object +#Example 1: using CST_RainFARM for a CSTools object nf <- 8 # Choose a downscaling by factor 8 exp <- 1 : (2 * 3 * 4 * 8 * 8) dim(exp) <- c(dataset = 1, member = 2, sdate = 3, ftime = 4, lat = 8, lon = 8) @@ -85,12 +88,25 @@ ww <- array(1., dim = c(8 * nf, 8 * nf)) res <- CST_RainFARM(data, nf, ww, nens=3) str(res) #List of 3 -# $ data: num [1:4, 1:64, 1:64, 1, 1:2, 1:3] 201 115 119 307 146 ... +# $ data: num [1, 1:2, 1:3, 1:3, 1:4, 1:64, 1:64] 260 553 281 278 143 ... # $ lon : num [1:64] 9.78 9.84 9.91 9.97 10.03 ... # $ lat : num [1:64] 39.8 39.8 39.9 40 40 ... dim(res$data) -# dataset member sdate ftime lat lon realization -# 1 2 3 4 64 64 3 +# dataset member realization sdate ftime lat lon +# 1 2 3 3 4 64 64 + +#Example 2: using CST_RainFARM for a CSTools object with parallel processing, +#dropping the "realization" dimension to be able to save results using +#\\code{CST_SaveExp}. +#Load dataset included in CSTools pacakge +data <- lonlat_prec +nf <- 8 +# Create a test array of weights +ww <- array(1., dim = c(dim(data$lon) * nf, dim(data$lat) * nf)) +res <- CST_RainFARM(data, nf, weights=ww, nprocs=2, drop_realization_dim = TRUE) +dim(res$data) +# dataset member sdate ftime lat lon +# 1 6 3 31 32 32 } \author{ Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} diff --git a/man/RainFARM.Rd b/man/RainFARM.Rd index 3ef2a2f6a4459d41b29d7d174763c3d4c8958404..984dcd42a6b7f374e3067ffe878ed19eff5fd181 100644 --- a/man/RainFARM.Rd +++ b/man/RainFARM.Rd @@ -5,8 +5,9 @@ \title{RainFARM stochastic precipitation downscaling (reduced version)} \usage{ RainFARM(data, lon, lat, nf, weights = 1, nens = 1, slope = 0, kmin = 1, - fglob = FALSE, fsmooth = TRUE, time_dim = NULL, lon_dim = "lon", - lat_dim = "lat", drop_realization_dim = FALSE, verbose = FALSE) + fglob = FALSE, fsmooth = TRUE, nprocs = 1, time_dim = NULL, + lon_dim = "lon", lat_dim = "lat", drop_realization_dim = FALSE, + verbose = FALSE) } \arguments{ \item{data}{Precipitation array to downscale. @@ -39,6 +40,8 @@ meaning that the slope is determined automatically over the dimensions specified \item{fsmooth}{Logical to conserve precipitation with a smoothing kernel (default: TRUE)} +\item{nprocs}{The number of parallel processes to spawn for the use for parallel computation in multiple cores. (default: 1)} + \item{time_dim}{String or character array with name(s) of time dimension(s) (e.g. "ftime", "sdate", "time" ...) over which to compute spectral slopes. If a character array of dimension names is provided, the spectral slopes @@ -51,6 +54,7 @@ is searched and the first one with more than one element is chosen.} \item{lat_dim}{Name of lat dimension ("lat" by default).} \item{drop_realization_dim}{Logical to remove the "realization" stochastic ensemble dimension (default: FALSE) + with the following behaviour if set to TRUE: 1) if \code{nens==1}: the dimension is dropped; diff --git a/tests/testthat/test-CST_RainFARM.R b/tests/testthat/test-CST_RainFARM.R new file mode 100644 index 0000000000000000000000000000000000000000..3608d4e8c3cd3347243289dd4b5d89cc4a9d0c41 --- /dev/null +++ b/tests/testthat/test-CST_RainFARM.R @@ -0,0 +1,122 @@ +context("Generic tests") +test_that("Sanity checks and simple use cases", { + # Generate simple synthetic data + # 4x5 in space, 2 members, 3 sdates, 6 ftime + r <- exp(rnorm(1 * 2 * 3 * 6 * 5 * 4)) + dim(r) <- c(dataset = 1, member = 2, sdate = 3, ftime = 6, lat = 5, lon = 4) + lon <- seq(0, 6, 2) + lat <- seq(10, 18, 2) + exp <- list(data = r, lat = lat, lon = lon) + attr(exp, 'class') <- 's2dv_cube' + + expect_warning( + res <- CST_RainFARM(exp, nf=8), + "input data are expected to be on a square grid with an even number of pixels per side." + ) + # Using a correct square dataset + r <- exp(rnorm(1 * 2 * 3 * 6 * 4 * 4)) + dim(r) <- c(dataset = 1, member = 2, sdate = 3, ftime = 6, lat = 4, lon = 4) + lat <- seq(10, 16, 2) + exp <- list(data = r, lat = lat, lon = lon) + attr(exp, 'class') <- 's2dv_cube' + + expect_warning( + res <- CST_RainFARM(exp, nf=8), + "Selected time dim: ftime" + ) + expect_error( + res <- CST_RainFARM(exp, nf=8, weights=array(0,dim=c(2,2))), + "The dimensions of the weights matrix" + ) + + dimexp=dim(exp$data) + + res <- CST_RainFARM(exp, nf=8, time_dim=c("ftime", "sdate"), slope=1.7, nens=2) + expect_equal(dim(res$data), c(dimexp["dataset"], dimexp["member"], + realization = 2, dimexp["sdate"], + dimexp["ftime"], dimexp["lat"] * 8, + dimexp["lon"] * 8)) + + expect_equivalent(length(res$lon), dimexp["lon"] * 8) + expect_equivalent(length(res$lat), dimexp["lat"] * 8) + + res <- CST_RainFARM(exp, nf=8, time_dim=c("ftime", "sdate"), + nens=2, drop_realization_dim=TRUE) + expect_equal(dim(res$data), c(dimexp["dataset"], dimexp["member"] * 2, + dimexp["sdate"], + dimexp["ftime"], dimexp["lat"] * 8, + dimexp["lon"] * 8)) + + res <- CST_RainFARM(exp, nf=8, time_dim=c("ftime", "sdate"), slope=1.7, + nens=2, nproc=2, fsmooth=FALSE) + expect_equal(dim(res$data), c(dimexp["dataset"], dimexp["member"], + realization = 2, dimexp["sdate"], + dimexp["ftime"], dimexp["lat"] * 8, + dimexp["lon"] * 8)) + + # Use rainfarmr package functions to test for consistency + + # Aggregated downscaled data should be the same as original + + expect_equivalent(agg(res$data[1,1,1,1,1,,], 4), exp$data[1,1,1,1,,]) + + res <- CST_RainFARM(exp, nf=8, time_dim=c("ftime", "sdate"), + nens=2, nproc=2, fglob=TRUE) + + expect_equal(mean(agg(res$data[1,1,1,1,1,,], 4)), + mean(exp$data[1,1,1,1,,])) + + # Create a more realistic perfect-model precipitation + z <- 1 : (32 * 32) + dim(z) <- c(32, 32) + x=rep(-15:16,32); dim(x) <- c(32, 32); x=x+4 ; y<-t(x); + z<-sqrt(x^2+y^2); z=(z<=12)*(z>=8)*0.8+0.6 + + exppf <- exp + nens = 10 + exppf$data <- array(1. , c(1, nens, 3, 6, 4, 4)) + dim(exppf$data) <- c(dataset = 1, member = nens, sdate = 3, + ftime = 6, lat = 4, lon = 4) + rpf <- array(1. , c(1, nens, 3, 6, 32, 32)) + for (i in 1:nens) { + for (j in 1:3) { + for (k in 1:6) { + f = initmetagauss(1.7, 32) + r = exp(metagauss(f)) * z + rpf[1, i, j, k, , ] <- r + exppf$data[1, i, j, k, , ] <- agg(r, 4) + } + } + } + rpfm=agg(apply(rpf, c(5, 6), mean),32) + + # Use climatological mean of PF precipitation to generate sythetic weights + w <- rfweights(rpfm, res$lon, res$lat, exp$lon, exp$lat, 8, fsmooth=FALSE ) + + res <- CST_RainFARM(exppf, nf=8, time_dim=c("ftime", "sdate", "member"), + nens=2, nproc=2, fsmooth=FALSE) + resw <- CST_RainFARM(exppf, nf=8, time_dim=c("ftime", "sdate", "member"), + nens=2, nproc=2, fsmooth=FALSE, weights=w) + resm <- agg(apply(res$data, c(6,7), mean),32) + reswm <- agg(apply(resw$data, c(6,7), mean),32) + + # Things should work much better with weights compared to no weights + expect_lt(mean((reswm-rpfm)^2), mean((resm-rpfm)^2)*0.5) + + # Test if downscaled data have similar slope to input data + # Synthetic fine-scale data with slope 1.7 + f = initmetagauss(1.7, 256) + r = exp(metagauss(f)) + ra <- agg(r, 8) + dim(ra) <- c(dataset = 1, member = 1, sdate = 1, ftime = 1, lat = 8, lon = 8) + expcoarse <- exp + expcoarse$data <- ra + dim(r) <- c(dataset = 1, member = 1, sdate = 1, ftime = 1, lat = 256, lon = 256) + expfine <- exp + expfine$data <- r + res <- CST_RainFARM(expcoarse, nf=32, time_dim=c("ftime", "sdate"), + slope=1.7, fsmooth=FALSE, drop_realization_dim=TRUE) + sres= CST_RFSlope(res, time_dim = c("ftime", "sdate")) + sexp= CST_RFSlope(expfine, time_dim = c("ftime", "sdate")) + expect_equal(sres, sexp, tolerance=0.25) +})