diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 48c1532058a450213bd230f8d7047211c986adc4..ee44142da4d9c2e6d181b9564e1d004263b83967 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -3,7 +3,8 @@ stages: build: stage: build script: - - module load R + - module load R/3.4.3-foss-2015a-bare Julia + - julia -e 'using Pkg; Pkg.add(["RainFARM"])' - R CMD build --resave-data . - - R CMD check --as-cran CSTools_*.tar.gz + - R CMD check --as-cran --no-manual CSTools_*.tar.gz - R -e 'covr::package_coverage()' diff --git a/DESCRIPTION b/DESCRIPTION index 781764ecfec0bfe0ad120a00f7b5d536b08e8aee..a579abe5e45e08b80081b6b443ee6a8e7157f86c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -4,17 +4,16 @@ Title: Set of Tools to Assess Skills Climate Forecasts on Seasonal-to-Decadal Version: 0.0.1 Authors@R: c( person("BSC-CNS", role = c("aut", "cph")), - person("Other copy rights holders", role = c("cph")), - person("Lauriane", "Batte", , "lauriane.batte@meteo.fr", role = "ctb"), person("Louis-Philippe", "Caron", , "louis-philippe.caron@bsc.es", role = "aut"), - person("Jost", "von Hardenberg", , "j.vonhardenberg@isac.cnr.it", role = "aut"), + person("Jost", "von Hardenberg", , "j.vonhardenberg@isac.cnr.it", role = c("aut", "cph")), + person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut", "cre")), person("Llorenç", "LLedo", , "llledo@bsc.es", role = "aut"), - person("Nicolau", "Manubens", , "nicolau.manubens@bsc.es", role = "ctb"), + person("Nicolau", "Manubens", , "nicolau.manubens@bsc.es", role = "aut"), person("Niti", "Mishra", , "niti.mishra@bsc.es", role = "aut"), - person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut", "cre")), - person("Bert", "Van Schaeybroeck", , "bertvs@meteo.be", role = "ctb"), person("Veronica", "Torralba", , "veronica.torralba@bsc.es", role = "aut"), - person("Deborah", "Verfaillie", , "deborah.verfaillie@bsc.es", role = "aut")) + person("Deborah", "Verfaillie", , "deborah.verfaillie@bsc.es", role = "aut"), + person("Lauriane", "Batte", , "lauriane.batte@meteo.fr", role = "ctb"), + person("Bert", "van Schaeybroeck", , "bertvs@meteo.be", role = "ctb")) Description: Improved global climate model calibration and regionalization techniques as well as better forecast verification methods need to be developed for Mediterranean region to extract as much climate information as possible from @@ -26,6 +25,8 @@ Imports: s2dverification, abind, stats, + JuliaCall, + ncdf4, ggplot2, plyr, data.table, @@ -34,8 +35,6 @@ Imports: Suggests: testthat License: Apache License 2.0 | file LICENSE -URL: https://earth.bsc.es/gitlab/external/cstools -BugReports: https://earth.bsc.es/gitlab/external/cstools/issues Encoding: UTF-8 LazyData: true -RoxygenNote: 6.0.1.9000 +RoxygenNote: 6.1.1 diff --git a/NAMESPACE b/NAMESPACE index 286f913e16af9efec932fbc5ff7b265244cc778f..dad92f85818b81a230118c79df6a4233c96702c6 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,10 +4,17 @@ export(CST_BiasCorrection) export(CST_Calibration) export(CST_MultiMetric) export(CST_MultivarRMSE) +export(CST_RFSlope) +export(CST_RainFARM) export(PlotForecastPDF) export(PlotMostLikelyQuantileMap) +export(RFSlope) +export(RFWeights) +export(RainFARM) +import(JuliaCall) import(ggplot2) import(multiApply) +import(ncdf4) import(s2dverification) import(stats) importFrom(data.table,CJ) diff --git a/R/CST_RFSlope.R b/R/CST_RFSlope.R new file mode 100644 index 0000000000000000000000000000000000000000..c1e10db42eb1f70ba3ab98ff0c04bf8fbc006d60 --- /dev/null +++ b/R/CST_RFSlope.R @@ -0,0 +1,164 @@ +#' @rdname CST_RFSlope +#' @title RainFARM spectral slopes from a CSTools object +#' +#' @author Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} +#' +#' @description This function computes spatial spectral slopes from a CSTools object +#' to be used for RainFARM stochastic precipitation downscaling method and accepts a CSTools object in input. +#' Please notice: The function uses the \code{julia_setup()} function from the JuliaCall library, +#' which takes a long time only the first time it is called on a machine. +#' +#' @param data CSTools object containing the spatial precipitation fields to downscale. +#' The data object is expected to have an element named \code{exp} with at least two +#' spatial dimensions named "lon" and "lat" and one or more dimensions over which +#' to average these slopes, which can be specified by parameter \code{time_dim}. +#' @param kmin First wavenumber for spectral slope (default \code{kmin=1}). +#' @param 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 +#' will be computed as an average over all elements belonging to those dimensions. +#' If omitted one of c("ftime", "sdate", "time") is searched and the first one with more than one element is chosen. +#' @return CST_RFSlope() returns spectral slopes using the RainFARM convention +#' (the logarithmic slope of k*|A(k)|^2 where A(k) are the spectral amplitudes). +#' The returned array has the same dimensions as the \code{exp} element of the input object, +#' minus the dimensions specified by \code{lon_dim}, \code{lat_dim} and \code{time_dim}. +#' @export +#' @examples +#' +#' #Example using CST_RFSlope for a CSTools object +#' exp <- 1 : (2 * 3 * 4 * 8 * 8) +#' dim(exp) <- c(dataset = 1, member = 2, sdate = 3, ftime = 4, lat = 8, lon = 8) +#' lon <- seq(10, 13.5, 0.5) +#' dim(lon) <- c(lon = length(lon)) +#' lat <- seq(40, 43.5, 0.5) +#' dim(lat) <- c(lat = length(lat)) +#' data <- list(exp = exp, lon = lon, lat = lat) +#' slopes <- CST_RFSlope(data) +#' dim(slopes) +#' # dataset member sdate +#' # 1 2 3 +#' slopes +#' # [,1] [,2] [,3] +#' #[1,] 1.893503 1.893503 1.893503 +#' #[2,] 1.893503 1.893503 1.893503 + +CST_RFSlope <- function(data, kmin = 1, time_dim = NULL) { + + slopes <- RFSlope(data$exp, kmin, time_dim, + lon_dim = "lon", lat_dim = "lat") + + return(slopes) +} + +#' @rdname RFSlope +#' @title RainFARM spectral slopes from an array (reduced version) +#' +#' @author Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} +#' +#' @description This function computes spatial spectral slopes from an array, +#' to be used for RainFARM stochastic precipitation downscaling method. +#' Please notice: The function uses the \code{julia_setup()} function from the JuliaCall library, +#' which takes a long time only the first time it is called on a machine. +#' +#' @param data Array containing the spatial precipitation fields to downscale. +#' The input array is expected to have at least two dimensions named "lon" and "lat" by default +#' (these default names can be changed with the \code{lon_dim} and \code{lat_dim} parameters) +#' and one or more dimensions over which to average the slopes, +#' which can be specified by parameter \code{time_dim}. +#' @param kmin First wavenumber for spectral slope (default \code{kmin=1}). +#' @param 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 +#' will be computed as an average over all elements belonging to those dimensions. +#' If omitted one of c("ftime", "sdate", "time") is searched and the first one +#' with more than one element is chosen. +#' @param lon_dim Name of lon dimension ("lon" by default). +#' @param lat_dim Name of lat dimension ("lat" by default). +#' @return RFSlope() returns spectral slopes using the RainFARM convention +#' (the logarithmic slope of k*|A(k)|^2 where A(k) are the spectral amplitudes). +#' The returned array has the same dimensions as the input array, +#' minus the dimensions specified by \code{lon_dim}, \code{lat_dim} and \code{time_dim}. +#' @export +#' @examples +#' # Example for the 'reduced' RFSlope function +#' +#' # Create a test array with dimension 8x8 and 20 timesteps, +#' # 3 starting dates and 20 ensemble members. +#' pr <- 1:(4*3*8*8*20) +#' dim(pr) <- c(ensemble = 4, sdate = 3, lon = 8, lat = 8, ftime = 20) +#' +#' # Compute the spectral slopes ignoring the wavenumber +#' # corresponding to the largest scale (the box) +#' slopes <- RFSlope(pr, kmin=2) +#' dim(slopes) +#' # ensemble sdate +#' # 4 3 +#' slopes +#' # [,1] [,2] [,3] +#' #[1,] 1.893503 1.893503 1.893503 +#' #[2,] 1.893503 1.893503 1.893503 +#' #[3,] 1.893503 1.893503 1.893503 +#' #[4,] 1.893503 1.893503 1.893503 + +RFSlope <- function(data, kmin = 1, time_dim = NULL, + lon_dim = "lon", lat_dim = "lat") { + + # Check/detect time_dim + if (is.null(time_dim)) { + time_dim_names <- c("ftime", "sdate", "time") + time_dim_num <- which(time_dim_names %in% names(dim(data))) + if (length(time_dim_num) > 0) { + # Find time dimension with length > 1 + ilong <- which(dim(data)[time_dim_names[time_dim_num]] > 1) + if (length(ilong) > 0) { + time_dim <- time_dim_names[time_dim_num[ilong[1]]] + } else { + stop("No time dimension longer than one found.") + } + } else { + stop("Could not automatically detect a target time dimension ", + "in the provided data in 'data'.") + } + print(paste("Selected time dim: ", time_dim)) + } + + # Check availability of JuliaCall + if (!("JuliaCall" %in% rownames(installed.packages()))) { + stop("Package 'JuliaCall' needs to be installed in order ", + "to run RainFARM.") + } + + # Perform common calls + JuliaCall::julia_setup(useRCall = FALSE) + julia_library("RainFARM") + + # reorder and group time_dim together at the end + cdim0 <- dim(data) + imask <- names(cdim0) %in% time_dim + data <- aperm(data,c(which(!imask), which(imask))) + cdim <- dim(data) + ind <- 1:length(which(!imask)) + # compact (multiply) time_dim dimensions + dim(data) <- c(cdim[ind], rainfarm_samples = prod(cdim[-ind])) + + # Repeatedly apply .RFSlope + print(dim(data)) + result <- Apply(data, c(lon_dim, lat_dim, "rainfarm_samples"), + .RFSlope, kmin)$output1 + + return(slopes = result) +} + +#' Atomic RFSlope +#' @param pr precipitation array to downscale with dims (lon, lat, time). +#' @param kmin first wavenumber for spectral slope (default kmin=1). +#' @return .RFSlope returns a scalar spectral slope using the RainFARM convention +#' (the logarithmic slope of k*|A(k)|^2 where A(k) is the spectral amplitude). +#' @noRd + +.RFSlope <- function(pr, kmin) { + + fxp <- julia_call("fft3d", pr, need_return = "R")[[1]] + sx <- julia_call("fitslopex", fxp, kmin = as.integer(kmin), + need_return = "R") + return(sx) +} diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R new file mode 100644 index 0000000000000000000000000000000000000000..144cd67fdb5914a095317912487f076ed007fee3 --- /dev/null +++ b/R/CST_RainFARM.R @@ -0,0 +1,298 @@ +#' @rdname CST_RainFARM +#' @title RainFARM stochastic precipitation downscaling of a CSTools object +#' +#' @author Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} +#' +#' @description This function implements the RainFARM stochastic precipitation +#' downscaling method and accepts a CSTools object in input. +#' Adapted for climate downscaling and including orographic correction +#' as described in Terzago et al. 2018. +#' References: Terzago, S. et al. (2018). NHESS 18(11), 2825–2840. +#' http://doi.org/10.5194/nhess-18-2825-2018 , +#' D'Onofrio et al. (2014), J of Hydrometeorology 15, 830-843; Rebora et. al. (2006), JHM 7, 724. +#' Please notice: The function uses the \code{julia_setup()} function from the JuliaCall library, +#' which takes a long time only the first time it is called on a machine. +#' +#' @param data CSTools object containing the spatial precipitation fields to downscale. +#' The data object is expected to have an element named \code{exp} with at least two +#' spatial dimensions named "lon" and "lat" and one or more dimensions over which +#' to compute average spectral slopes (unless specified with parameter \code{slope}), +#' which can be specified by parameter \code{time_dim}. +#' @param weights Matrix with climatological weights which can be obtained using +#' the \code{RFWeights} function. If \code{weights=1.} (default) no weights are used. +#' The matrix should have dimensions (lon, lat) in this order. +#' The names of these dimensions are not checked. +#' @param nf Refinement factor for downscaling (the output resolution is increased by this factor). +#' @param slope Prescribed spectral slope. The default is \code{slope=0.} +#' meaning that the slope is determined automatically over the dimensions specified by \code{time_dim}. +#' @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 conserve global precipitation over the domain (default: FALSE). +#' @param fsmooth Logical to conserve precipitation with a smoothing kernel (default: TRUE). +#' @param 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 +#' will be computed as an average over all elements belonging to those dimensions. +#' 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 of the Julia package (default: FALSE). +#' @param 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; +#' +#' 2) if \code{nens>1} and a "member" dimension exists: +#' 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". +#' +#' @return CST_RainFARM() returns a downscaled CSTools object. +#' If \code{nens>1} an additional dimension named "realization" is added to the \code{exp} array +#' after the "member" dimension (unless \code{drop_realization_dim=TRUE} is specified). +#' The ordering of the remaining dimensions in the \code{exp} element of the input object is maintained. +#' @import multiApply +#' @export +#' @examples +#' #Example 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) +#' lon <- seq(10, 13.5, 0.5) +#' dim(lon) <- c(lon = length(lon)) +#' lat <- seq(40, 43.5, 0.5) +#' dim(lat) <- c(lat = length(lat)) +#' data <- list(exp = exp, lon = lon, lat = lat) +#' # Create a test array of weights +#' ww <- array(1., dim = c(8 * nf, 8 * nf)) +#' res <- CST_RainFARM(data, nf, ww, nens=3) +#' str(res) +#' #List of 3 +#' # $ exp: num [1:4, 1:64, 1:64, 1, 1:2, 1:3] 201 115 119 307 146 ... +#' # $ 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$exp) +#' # dataset member sdate ftime lat lon realization +#' # 1 2 3 4 64 64 3 +#' +CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1, + nens = 1, fglob = FALSE, fsmooth = TRUE, + time_dim = NULL, verbose = FALSE, + drop_realization_dim = FALSE) { + + res <- RainFARM(data$exp, data$lon, data$lat, + nf, weights, nens, slope, kmin, fglob, fsmooth, + time_dim, lon_dim = "lon", lat_dim = "lat", + drop_realization_dim, verbose) + + data$exp <- res$data + data$lon <- res$lon + data$lat <- res$lat + + return(data) +} + +#' @rdname RainFARM +#' @title RainFARM stochastic precipitation downscaling (reduced version) +#' @author Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} +#' @description This function implements the RainFARM stochastic precipitation downscaling method +#' and accepts in input an array with named dims ("lon", "lat") +#' and one or more dimension (such as "ftime", "sdate" or "time") +#' over which to average automatically determined spectral slopes. +#' Adapted for climate downscaling and including orographic correction. +#' References: +#' Terzago, S. et al. (2018). NHESS 18(11), 2825–2840. http://doi.org/10.5194/nhess-18-2825-2018, +#' D'Onofrio et al. (2014), J of Hydrometeorology 15, 830-843; Rebora et. al. (2006), JHM 7, 724. +#' Please notice: The function uses the \code{julia_setup()} function from the JuliaCall library, +#' which takes a long time only the first time it is called on a machine. +#' @param data Precipitation array to downscale. +#' The input array is expected to have at least two dimensions named "lon" and "lat" by default +#' (these default names can be changed with the \code{lon_dim} and \code{lat_dim} parameters) +#' and one or more dimensions over which to average these slopes, +#' which can be specified by parameter \code{time_dim}. +#' @param lon Vector or array of longitudes. +#' @param lat Vector or array of latitudes. +#' @param weights Matrix with climatological weights which can be obtained using +#' the \code{RFWeights} function. If \code{weights=1.} (default) no weights are used. +#' The matrix should have dimensions (lon, lat) in this order. +#' The names of these dimensions are not checked. +#' @param nf Refinement factor for downscaling (the output resolution is increased by this factor). +#' @param slope Prescribed spectral slope. The default is \code{slope=0.} +#' meaning that the slope is determined automatically over the dimensions specified by \code{time_dim}. +#' @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 smoothing kernel (default: TRUE) +#' @param 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 +#' will be computed over all elements belonging to those dimensions. +#' If omitted one of c("ftime", "sdate", "time") +#' is searched and the first one with more than one element is chosen. +#' @param lon_dim Name of lon dimension ("lon" by default). +#' @param lat_dim Name of lat dimension ("lat" by default). +#' @param verbose logical for verbose output of the Julia package (default: FALSE). +#' @param 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; +#' +#' 2) if \code{nens>1} and a "member" dimension exists: +#' 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". +#' +#' @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 +#' after the "member" dimension (if it exists and unless \code{drop_realization_dim=TRUE} is specified). +#' The ordering of the remaining dimensions in the \code{exp} element of the input object is maintained. +#' @import JuliaCall +#' @import multiApply +#' @export +#' @examples +#' # Example for the 'reduced' RainFARM function +#' nf <- 8 # Choose a downscaling by factor 8 +#' nens <- 3 # Number of ensemble members +#' # create a test array with dimension 8x8 and 20 timesteps +#' # or provide your own read from a netcdf file +#' pr <- rnorm(8 * 8 * 20) +#' dim(pr) <- c(lon = 8, lat = 8, ftime = 20) +#' lon_mat <- seq(10, 13.5, 0.5) # could also be a 2d matrix +#' lat_mat <- seq(40, 43.5, 0.5) +#' # Create a test array of weights +#' ww <- array(1., dim = c(8 * nf, 8 * nf)) +#' # ... or create proper weights using an external fine-scale climatology file +#' # Specify a weightsfn filename if you wish to save the weights +#' \donttest{ +#' ww <- RFWeights("./worldclim.nc", nf, lon = lon_mat, lat = lat_mat, +#' weightsfn = "", fsmooth = TRUE) +#' } +#' # downscale using weights (ww=1. means do not use weights) +#' res <- RainFARM(pr, lon_mat, lat_mat, nf, +#' fsmooth = TRUE, fglob = FALSE, +#' weights = ww, nens = 2, verbose = TRUE) +#' str(res) +#' #List of 3 +#' # $ data: num [1:3, 1:20, 1:64, 1:64] 0.186 0.212 0.138 3.748 0.679 ... +#' # $ 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) +#' # lon lat ftime realization +#' # 64 64 20 2 +#' +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", + drop_realization_dim = FALSE, verbose = FALSE) { + + # Check/detect time_dim + if (is.null(time_dim)) { + time_dim_names <- c("ftime", "sdate", "time") + time_dim_num <- which(time_dim_names %in% names(dim(data))) + if (length(time_dim_num) > 0) { + # Find time dimension with length > 1 + ilong <- which(dim(data)[time_dim_names[time_dim_num]] > 1) + if (length(ilong) > 0) { + time_dim <- time_dim_names[time_dim_num[ilong[1]]] + } else { + stop("No time dimension longer than one found.") + } + } else { + stop("Could not automatically detect a target time dimension ", + "in the provided data in 'data'.") + } + print(paste("Selected time dim: ", time_dim)) + } + + # Check availability of JuliaCall + if (!("JuliaCall" %in% rownames(installed.packages()))) { + stop("Package 'JuliaCall' needs to be installed in order ", + "to run RainFARM.") + } + + # Perform common calls + JuliaCall::julia_setup(useRCall = FALSE) + julia_library("RainFARM") + + r <- julia_call("lon_lat_fine", lon, lat, as.integer(nf), + need_return = "R") + lon_f <- r[[1]] + lat_f <- r[[2]] + + # reorder and group time_dim together at the end + cdim0 <- dim(data) + imask <- names(cdim0) %in% time_dim + data <- aperm(data, c(which(!imask), which(imask))) + cdim <- dim(data) + ind <- 1:length(which(!imask)) + # compact (multiply) time_dim dimensions + dim(data) <- c(cdim[ind], rainfarm_samples = prod(cdim[-ind])) + + # Repeatedly apply .RainFARM + result <- Apply(data, c(lon_dim, lat_dim, "rainfarm_samples"), .RainFARM, + weights, nf, nens, slope, kmin, + fglob, fsmooth, verbose)$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)]) + + # 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)) + iorder <- sapply(c(names(cdim0)[1:ienspos], "realization", names(cdim0)[-(1:ienspos)]), grep, names(dim(result))) + result <- aperm(result, iorder) + + if (drop_realization_dim) { + cdim = dim(result) + if (nens == 1) { + # just drop it if only one member + dim(result) <- cdim[-which(names(cdim) == "realization")[1]] + } else if("member" %in% names(cdim)) { + # compact member and realization dimension if member dim exists, else rename realization to member + ind <- which(names(cdim) %in% c("member", "realization")) + dim(result) <- c(cdim[1:(ind[1]-1)], cdim[ind[1]]*cdim[ind[2]], cdim[(ind[2]+1):length(cdim)]) + } else { + ind <- which(names(cdim) %in% "realization") + names(dim(result))[ind] <- "member" + } + } + + return(list(data = result, lon = lon_f, lat = lat_f)) +} + +#' Atomic RainFARM +#' @param pr Precipitation array to downscale with dimensions (lon, lat, time). +#' @param weights Matrix with climatological weights which can be obtained using +#' the \code{RFWeights} function (default: \code{weights=1.} i.e. no weights). +#' @param nf Refinement factor for downscaling (the output resolution is increased by this factor). +#' @param slope Prescribed spectral slope (default: \code{slope=0.} +#' meaning that the slope is determined automatically over the dimensions specified by \code{time_dim}. +#' @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 verbose Logical for verbose output of the Julia package (default: FALSE). +#' @return .RainFARM returns a downscaled array with dimensions (lon, lat, time, realization) +#' @noRd +.RainFARM <- function(pr, weights, nf, nens, slope, kmin, + fglob, fsmooth, verbose) { + + if (slope == 0) { + fxp <- julia_call("fft3d", pr, need_return = "R")[[1]] + sx <- julia_call("fitslopex", fxp, kmin = as.integer(kmin), + need_return = "R") + } else { + sx <- slope + } + + result_dims <- c(dim(pr)[1] * nf, dim(pr)[2] * nf, dim(pr)[3], + realization = nens) + r <- array(dim = result_dims) + + for (i in 1:nens) { + r[, , , i] <- julia_call("rainfarm", pr, sx, as.integer(nf), weights, + fglob = fglob, fsmooth = fsmooth, + verbose = verbose, need_return = "R") + } + return(r) +} diff --git a/R/RFWeights.R b/R/RFWeights.R new file mode 100644 index 0000000000000000000000000000000000000000..66b592e187893a6f0d2f63805b8138f6c7dd345c --- /dev/null +++ b/R/RFWeights.R @@ -0,0 +1,77 @@ +#' Compute climatological weights for RainFARM stochastic precipitation downscaling +#' +#' @author Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} +#' +#' @description Compute climatological ("orographic") weights from a fine-scale precipitation climatology file. +#' Reference: Terzago, S., Palazzi, E., & von Hardenberg, J. (2018). +#' Stochastic downscaling of precipitation in complex orography: +#' A simple method to reproduce a realistic fine-scale climatology. +#' Natural Hazards and Earth System Sciences, 18(11), +#' 2825–2840. http://doi.org/10.5194/nhess-18-2825-2018 . +#' Please notice: The function uses the \code{julia_setup()} +#' function from the JuliaCall library, +#' which takes a long time only the first time it is called on a machine +#' @param orofile Filename of a fine-scale precipitation climatology. +#' The file is expected to be in NetCDF format and should contain +#' at least one precipitation field. If several fields at different times are provided, +#' a climatology is derived by time averaging. +#' Suitable climatology files could be for example a fine-scale precipitation climatology +#' from a high-resolution regional climate model (see e.g. Terzago et al. 2018), a local +#' high-resolution gridded climatology from observations, or a reconstruction such as those which +#' can be downloaded from the WORLDCLIM (http://www.worldclim.org) or CHELSA (http://chelsa-climate.org) +#' websites. The latter data will need to be converted to NetCDF format before being used (see for example the GDAL tools (https://www.gdal.org). +#' @param nf Refinement factor for downscaling (the output resolution is increased by this factor). +#' @param lon Vector or array of longitudes (can be omitted if \code{reffile} is used). +#' @param lat Vector or array of latitudes (can be omitted if \code{reffile} is used). +#' @param reffile Filename of reference file (e.g. the file to downscale). +#' Not needed if \code{lon} and \code{lat} are passed. +#' @param varname Name of the variable to be read from \code{reffile}. +#' @param weightsfn Optional filename to save the weights also to an external file in netcdf format. +#' @param fsmooth Logical to use smooth conservation (default) or large-scale box-average conservation. +#' @return A matrix containing the weights with dimensions (lon, lat). +#' @import ncdf4 +#' @import JuliaCall +#' @examples +#' # Create weights to be used with the CST_RainFARM() or RainFARM() functions +#' # using an external fine-scale climatology file. +#' +#' \donttest{ +#' # Specify lon and lat of the input +#' lon_mat <- seq(10,13.5,0.5) # could also be a 2d matrix +#' lat_mat <- seq(40,43.5,0.5) +#' +#' nf <- 8 +#' ww <- RFWeights("./worldclim.nc", nf, lon=lon_mat, lat=lat_mat, weightsfn="", fsmooth = TRUE) +#' +#' # or use a reference file (e.g. the input file itself) to specify the coordinates +#' ww <- RFWeights("./worldclim.nc", nf, reffile="./input.nc") +#' } +#' @export + +RFWeights <- function(orofile, nf, lon=NULL, lat=NULL, reffile="", + weightsfn="", varname="", fsmooth=TRUE) { + + julia_setup(useRCall = FALSE) + julia_library("RainFARM") + if ( (reffile == "") & ( (typeof(lon) == "NULL") | (typeof(lat) == "NULL"))) { + stop("Neither a reference filename nor lon and lat specified!\n") + } + if (reffile == "") { + xvar <- ncdim_def("lon", "degrees_east", lon, longname = "longitude") + yvar <- ncdim_def("lat", "degrees_north", lat, longname = "latitude") + tvar <- ncdim_def("time", "days since 01-01-2000", 0., longname = "time", + unlim = T) + my_ncdf <- ncvar_def("grid", "m", list(xvar, yvar, tvar), -999, + longname = "grid", prec = "single") + reffile <- tempfile() + ncfile <- nc_create(reffile, my_ncdf, force_v4 = FALSE) + nc_close(ncfile) + } + + ww <- julia_call("rfweights", orofile, reffile, as.integer(nf), + weightsfn = weightsfn, varname = varname, fsmooth = fsmooth) + + attributes(dim(ww))$names <- c("lon","lat") + unlink(reffile) + return(ww) +} diff --git a/man/CST_RFSlope.Rd b/man/CST_RFSlope.Rd new file mode 100644 index 0000000000000000000000000000000000000000..a931374068f7ce9b762f76641faae54605e4420d --- /dev/null +++ b/man/CST_RFSlope.Rd @@ -0,0 +1,55 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_RFSlope.R +\name{CST_RFSlope} +\alias{CST_RFSlope} +\title{RainFARM spectral slopes from a CSTools object} +\usage{ +CST_RFSlope(data, kmin = 1, time_dim = NULL) +} +\arguments{ +\item{data}{CSTools object containing the spatial precipitation fields to downscale. +The data object is expected to have an element named \code{exp} with at least two +spatial dimensions named "lon" and "lat" and one or more dimensions over which +to average these slopes, which can be specified by parameter \code{time_dim}.} + +\item{kmin}{First wavenumber for spectral slope (default \code{kmin=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 +will be computed as an average over all elements belonging to those dimensions. +If omitted one of c("ftime", "sdate", "time") is searched and the first one with more than one element is chosen.} +} +\value{ +CST_RFSlope() returns spectral slopes using the RainFARM convention +(the logarithmic slope of k*|A(k)|^2 where A(k) are the spectral amplitudes). +The returned array has the same dimensions as the \code{exp} element of the input object, +minus the dimensions specified by \code{lon_dim}, \code{lat_dim} and \code{time_dim}. +} +\description{ +This function computes spatial spectral slopes from a CSTools object +to be used for RainFARM stochastic precipitation downscaling method and accepts a CSTools object in input. +Please notice: The function uses the \code{julia_setup()} function from the JuliaCall library, +which takes a long time only the first time it is called on a machine. +} +\examples{ + +#Example using CST_RFSlope for a CSTools object +exp <- 1 : (2 * 3 * 4 * 8 * 8) +dim(exp) <- c(dataset = 1, member = 2, sdate = 3, ftime = 4, lat = 8, lon = 8) +lon <- seq(10, 13.5, 0.5) +dim(lon) <- c(lon = length(lon)) +lat <- seq(40, 43.5, 0.5) +dim(lat) <- c(lat = length(lat)) +data <- list(exp = exp, lon = lon, lat = lat) +slopes <- CST_RFSlope(data) +dim(slopes) +# dataset member sdate +# 1 2 3 +slopes +# [,1] [,2] [,3] +#[1,] 1.893503 1.893503 1.893503 +#[2,] 1.893503 1.893503 1.893503 +} +\author{ +Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} +} diff --git a/man/CST_RainFARM.Rd b/man/CST_RainFARM.Rd new file mode 100644 index 0000000000000000000000000000000000000000..2120194e1eee8ac44ef3525fbf01e9511b0cbd16 --- /dev/null +++ b/man/CST_RainFARM.Rd @@ -0,0 +1,97 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_RainFARM.R +\name{CST_RainFARM} +\alias{CST_RainFARM} +\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) +} +\arguments{ +\item{data}{CSTools object containing the spatial precipitation fields to downscale. +The data object is expected to have an element named \code{exp} with at least two +spatial dimensions named "lon" and "lat" and one or more dimensions over which +to compute average spectral slopes (unless specified with parameter \code{slope}), +which can be specified by parameter \code{time_dim}.} + +\item{nf}{Refinement factor for downscaling (the output resolution is increased by this factor).} + +\item{weights}{Matrix with climatological weights which can be obtained using +the \code{RFWeights} function. If \code{weights=1.} (default) no weights are used. +The matrix should have dimensions (lon, lat) in this order. +The names of these dimensions are not checked.} + +\item{slope}{Prescribed spectral slope. The default is \code{slope=0.} +meaning that the slope is determined automatically over the dimensions specified by \code{time_dim}.} + +\item{kmin}{First wavenumber for spectral slope (default: \code{kmin=1}).} + +\item{nens}{Number of ensemble members to produce (default: \code{nens=1}).} + +\item{fglob}{Logical to conserve global precipitation over the domain (default: FALSE).} + +\item{fsmooth}{Logical to conserve precipitation with a smoothing kernel (default: TRUE).} + +\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 +will be computed as an average over all elements belonging to those dimensions. +If omitted one of c("ftime", "sdate", "time") is searched and the first one with more +than one element is chosen.} + +\item{verbose}{Logical for verbose output of the Julia package (default: FALSE).} + +\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; + +2) if \code{nens>1} and a "member" dimension exists: + 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".} +} +\value{ +CST_RainFARM() returns a downscaled CSTools object. +If \code{nens>1} an additional dimension named "realization" is added to the \code{exp} array +after the "member" dimension (unless \code{drop_realization_dim=TRUE} is specified). +The ordering of the remaining dimensions in the \code{exp} element of the input object is maintained. +} +\description{ +This function implements the RainFARM stochastic precipitation +downscaling method and accepts a CSTools object in input. +Adapted for climate downscaling and including orographic correction +as described in Terzago et al. 2018. +References: Terzago, S. et al. (2018). NHESS 18(11), 2825–2840. +http://doi.org/10.5194/nhess-18-2825-2018 , +D'Onofrio et al. (2014), J of Hydrometeorology 15, 830-843; Rebora et. al. (2006), JHM 7, 724. +Please notice: The function uses the \code{julia_setup()} function from the JuliaCall library, +which takes a long time only the first time it is called on a machine. +} +\examples{ +#Example 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) +lon <- seq(10, 13.5, 0.5) +dim(lon) <- c(lon = length(lon)) +lat <- seq(40, 43.5, 0.5) +dim(lat) <- c(lat = length(lat)) +data <- list(exp = exp, lon = lon, lat = lat) +# Create a test array of weights +ww <- array(1., dim = c(8 * nf, 8 * nf)) +res <- CST_RainFARM(data, nf, ww, nens=3) +str(res) +#List of 3 +# $ exp: num [1:4, 1:64, 1:64, 1, 1:2, 1:3] 201 115 119 307 146 ... +# $ 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$exp) +# dataset member sdate ftime lat lon realization +# 1 2 3 4 64 64 3 + +} +\author{ +Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} +} diff --git a/man/PlotForecastPDF.Rd b/man/PlotForecastPDF.Rd index bd273aafe178ac7664fe4090def7fb752e6857bf..0af93c072220be143acb582921d3aff2a23a5330 100644 --- a/man/PlotForecastPDF.Rd +++ b/man/PlotForecastPDF.Rd @@ -4,10 +4,11 @@ \alias{PlotForecastPDF} \title{Plot one or multiple ensemble forecast pdfs for the same event} \usage{ -PlotForecastPDF(fcst, tercile.limits, extreme.limits = NULL, obs = NULL, - plotfile = NULL, title = "Set a title", var.name = "Varname (units)", - fcst.names = NULL, add.ensmemb = c("above", "below", "no"), - color.set = c("ggplot", "s2s4e", "hydro")) +PlotForecastPDF(fcst, tercile.limits, extreme.limits = NULL, + obs = NULL, plotfile = NULL, title = "Set a title", + var.name = "Varname (units)", fcst.names = NULL, + add.ensmemb = c("above", "below", "no"), color.set = c("ggplot", + "s2s4e", "hydro")) } \arguments{ \item{fcst}{a dataframe or array containing all the ensember members for each frecast. If \code{'fcst'} is an array, it should have two labelled dimensions, and one of them should be \code{'members'}. If \code{'fcsts'} is a data.frame, each column shoul be a separate forecast, with the rows beeing the different ensemble members.} diff --git a/man/RFSlope.Rd b/man/RFSlope.Rd new file mode 100644 index 0000000000000000000000000000000000000000..1f5689aac73e0a5b2489bf25eb0bf2cc98412e9d --- /dev/null +++ b/man/RFSlope.Rd @@ -0,0 +1,65 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_RFSlope.R +\name{RFSlope} +\alias{RFSlope} +\title{RainFARM spectral slopes from an array (reduced version)} +\usage{ +RFSlope(data, kmin = 1, time_dim = NULL, lon_dim = "lon", + lat_dim = "lat") +} +\arguments{ +\item{data}{Array containing the spatial precipitation fields to downscale. +The input array is expected to have at least two dimensions named "lon" and "lat" by default +(these default names can be changed with the \code{lon_dim} and \code{lat_dim} parameters) +and one or more dimensions over which to average the slopes, +which can be specified by parameter \code{time_dim}.} + +\item{kmin}{First wavenumber for spectral slope (default \code{kmin=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 +will be computed as an average over all elements belonging to those dimensions. +If omitted one of c("ftime", "sdate", "time") is searched and the first one +with more than one element is chosen.} + +\item{lon_dim}{Name of lon dimension ("lon" by default).} + +\item{lat_dim}{Name of lat dimension ("lat" by default).} +} +\value{ +RFSlope() returns spectral slopes using the RainFARM convention +(the logarithmic slope of k*|A(k)|^2 where A(k) are the spectral amplitudes). +The returned array has the same dimensions as the input array, +minus the dimensions specified by \code{lon_dim}, \code{lat_dim} and \code{time_dim}. +} +\description{ +This function computes spatial spectral slopes from an array, +to be used for RainFARM stochastic precipitation downscaling method. +Please notice: The function uses the \code{julia_setup()} function from the JuliaCall library, +which takes a long time only the first time it is called on a machine. +} +\examples{ +# Example for the 'reduced' RFSlope function + +# Create a test array with dimension 8x8 and 20 timesteps, +# 3 starting dates and 20 ensemble members. +pr <- 1:(4*3*8*8*20) +dim(pr) <- c(ensemble = 4, sdate = 3, lon = 8, lat = 8, ftime = 20) + +# Compute the spectral slopes ignoring the wavenumber +# corresponding to the largest scale (the box) +slopes <- RFSlope(pr, kmin=2) +dim(slopes) +# ensemble sdate +# 4 3 +slopes +# [,1] [,2] [,3] +#[1,] 1.893503 1.893503 1.893503 +#[2,] 1.893503 1.893503 1.893503 +#[3,] 1.893503 1.893503 1.893503 +#[4,] 1.893503 1.893503 1.893503 +} +\author{ +Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} +} diff --git a/man/RFWeights.Rd b/man/RFWeights.Rd new file mode 100644 index 0000000000000000000000000000000000000000..1f90924f2f5aaa0dd9197c11d4c5d01f799cf652 --- /dev/null +++ b/man/RFWeights.Rd @@ -0,0 +1,68 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RFWeights.R +\name{RFWeights} +\alias{RFWeights} +\title{Compute climatological weights for RainFARM stochastic precipitation downscaling} +\usage{ +RFWeights(orofile, nf, lon = NULL, lat = NULL, reffile = "", + weightsfn = "", varname = "", fsmooth = TRUE) +} +\arguments{ +\item{orofile}{Filename of a fine-scale precipitation climatology. +The file is expected to be in NetCDF format and should contain +at least one precipitation field. If several fields at different times are provided, +a climatology is derived by time averaging. +Suitable climatology files could be for example a fine-scale precipitation climatology +from a high-resolution regional climate model (see e.g. Terzago et al. 2018), a local +high-resolution gridded climatology from observations, or a reconstruction such as those which +can be downloaded from the WORLDCLIM (http://www.worldclim.org) or CHELSA (http://chelsa-climate.org) +websites. The latter data will need to be converted to NetCDF format before being used (see for example the GDAL tools (https://www.gdal.org).} + +\item{nf}{Refinement factor for downscaling (the output resolution is increased by this factor).} + +\item{lon}{Vector or array of longitudes (can be omitted if \code{reffile} is used).} + +\item{lat}{Vector or array of latitudes (can be omitted if \code{reffile} is used).} + +\item{reffile}{Filename of reference file (e.g. the file to downscale). +Not needed if \code{lon} and \code{lat} are passed.} + +\item{weightsfn}{Optional filename to save the weights also to an external file in netcdf format.} + +\item{varname}{Name of the variable to be read from \code{reffile}.} + +\item{fsmooth}{Logical to use smooth conservation (default) or large-scale box-average conservation.} +} +\value{ +A matrix containing the weights with dimensions (lon, lat). +} +\description{ +Compute climatological ("orographic") weights from a fine-scale precipitation climatology file. +Reference: Terzago, S., Palazzi, E., & von Hardenberg, J. (2018). +Stochastic downscaling of precipitation in complex orography: +A simple method to reproduce a realistic fine-scale climatology. +Natural Hazards and Earth System Sciences, 18(11), +2825–2840. http://doi.org/10.5194/nhess-18-2825-2018 . +Please notice: The function uses the \code{julia_setup()} +function from the JuliaCall library, +which takes a long time only the first time it is called on a machine +} +\examples{ +# Create weights to be used with the CST_RainFARM() or RainFARM() functions +# using an external fine-scale climatology file. + +\donttest{ +# Specify lon and lat of the input +lon_mat <- seq(10,13.5,0.5) # could also be a 2d matrix +lat_mat <- seq(40,43.5,0.5) + +nf <- 8 +ww <- RFWeights("./worldclim.nc", nf, lon=lon_mat, lat=lat_mat, weightsfn="", fsmooth = TRUE) + +# or use a reference file (e.g. the input file itself) to specify the coordinates +ww <- RFWeights("./worldclim.nc", nf, reffile="./input.nc") +} +} +\author{ +Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} +} diff --git a/man/RainFARM.Rd b/man/RainFARM.Rd new file mode 100644 index 0000000000000000000000000000000000000000..b02b4102fe3868513c41aaccf41ba0f70086d64b --- /dev/null +++ b/man/RainFARM.Rd @@ -0,0 +1,116 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_RainFARM.R +\name{RainFARM} +\alias{RainFARM} +\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) +} +\arguments{ +\item{data}{Precipitation array to downscale. +The input array is expected to have at least two dimensions named "lon" and "lat" by default +(these default names can be changed with the \code{lon_dim} and \code{lat_dim} parameters) +and one or more dimensions over which to average these slopes, +which can be specified by parameter \code{time_dim}.} + +\item{lon}{Vector or array of longitudes.} + +\item{lat}{Vector or array of latitudes.} + +\item{nf}{Refinement factor for downscaling (the output resolution is increased by this factor).} + +\item{weights}{Matrix with climatological weights which can be obtained using +the \code{RFWeights} function. If \code{weights=1.} (default) no weights are used. +The matrix should have dimensions (lon, lat) in this order. +The names of these dimensions are not checked.} + +\item{nens}{Number of ensemble members to produce (default: \code{nens=1}).} + +\item{slope}{Prescribed spectral slope. The default is \code{slope=0.} +meaning that the slope is determined automatically over the dimensions specified by \code{time_dim}.} + +\item{kmin}{First wavenumber for spectral slope (default: \code{kmin=1}).} + +\item{fglob}{Logical to conseve global precipitation over the domain (default: FALSE)} + +\item{fsmooth}{Logical to conserve precipitation with a smoothing kernel (default: TRUE)} + +\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 +will be computed over all elements belonging to those dimensions. +If omitted one of c("ftime", "sdate", "time") +is searched and the first one with more than one element is chosen.} + +\item{lon_dim}{Name of lon dimension ("lon" by default).} + +\item{lat_dim}{Name of lat dimension ("lat" by default).} + +\item{drop_realization_dim}{Logical to remove the "realization" stochastic ensemble dimension (default: FALSE) + +1) if \code{nens==1}: the dimension is dropped; + +2) if \code{nens>1} and a "member" dimension exists: + 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".} + +\item{verbose}{logical for verbose output of the Julia package (default: FALSE).} +} +\value{ +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 +after the "member" dimension (if it exists and unless \code{drop_realization_dim=TRUE} is specified). +The ordering of the remaining dimensions in the \code{exp} element of the input object is maintained. +} +\description{ +This function implements the RainFARM stochastic precipitation downscaling method +and accepts in input an array with named dims ("lon", "lat") +and one or more dimension (such as "ftime", "sdate" or "time") +over which to average automatically determined spectral slopes. +Adapted for climate downscaling and including orographic correction. +References: +Terzago, S. et al. (2018). NHESS 18(11), 2825–2840. http://doi.org/10.5194/nhess-18-2825-2018, +D'Onofrio et al. (2014), J of Hydrometeorology 15, 830-843; Rebora et. al. (2006), JHM 7, 724. +Please notice: The function uses the \code{julia_setup()} function from the JuliaCall library, +which takes a long time only the first time it is called on a machine. +} +\examples{ +# Example for the 'reduced' RainFARM function +nf <- 8 # Choose a downscaling by factor 8 +nens <- 3 # Number of ensemble members +# create a test array with dimension 8x8 and 20 timesteps +# or provide your own read from a netcdf file +pr <- rnorm(8 * 8 * 20) +dim(pr) <- c(lon = 8, lat = 8, ftime = 20) +lon_mat <- seq(10, 13.5, 0.5) # could also be a 2d matrix +lat_mat <- seq(40, 43.5, 0.5) +# Create a test array of weights +ww <- array(1., dim = c(8 * nf, 8 * nf)) +# ... or create proper weights using an external fine-scale climatology file +# Specify a weightsfn filename if you wish to save the weights +\donttest{ +ww <- RFWeights("./worldclim.nc", nf, lon = lon_mat, lat = lat_mat, + weightsfn = "", fsmooth = TRUE) +} +# downscale using weights (ww=1. means do not use weights) +res <- RainFARM(pr, lon_mat, lat_mat, nf, + fsmooth = TRUE, fglob = FALSE, + weights = ww, nens = 2, verbose = TRUE) +str(res) +#List of 3 +# $ data: num [1:3, 1:20, 1:64, 1:64] 0.186 0.212 0.138 3.748 0.679 ... +# $ 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) +# lon lat ftime realization +# 64 64 20 2 + +} +\author{ +Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} +}