diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index ee44142da4d9c2e6d181b9564e1d004263b83967..866209853288e1301b4b6787c5f70af2055b7249 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -3,8 +3,9 @@ stages: build: stage: build script: - - module load R/3.4.3-foss-2015a-bare Julia - - julia -e 'using Pkg; Pkg.add(["RainFARM"])' + #- module load R/3.4.3-foss-2015a-bare Julia + - module load R + #- julia -e 'using Pkg; Pkg.add(["RainFARM"])' - R CMD build --resave-data . - - R CMD check --as-cran --no-manual CSTools_*.tar.gz + - R CMD check --as-cran --no-manual --run-donttest CSTools_*.tar.gz - R -e 'covr::package_coverage()' diff --git a/DESCRIPTION b/DESCRIPTION index 6b64acde7b126e43fcf2358e2fcabe09cdb5f895..f4ad9b9fab9e2f60d72889f6f144952f571e3e05 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,10 +1,11 @@ Package: CSTools -Title: Set of Tools to Assess Skill of Climate Forecasts on Seasonal-to-Decadal Timescales +Title: Set of Tools to Assess Skill of Climate Forecasts on Seasonal-to-Decadal + Timescales Version: 1.0.0 Authors@R: c( person("BSC-CNS", role = c("cph")), person("Louis-Philippe", "Caron", , "louis-philippe.caron@bsc.es", role = "aut", comment = c(ORCID = "0000-0001-5221-0147")), - person("Jost", "von Hardenberg", , "j.vonhardenberg@isac.cnr.it", role = c("aut", "cph")), + person("Jost", "von Hardenberg", , "j.vonhardenberg@isac.cnr.it", role = "aut", comment = c(ORCID = "0000-0002-5312-8070")), 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 = "aut"), @@ -14,24 +15,29 @@ Authors@R: c( person("Lauriane", "Batte", , "lauriane.batte@meteo.fr", role = "ctb"), person("Jesus", "Peña", , "jesus.pena@bsc.es", role = "ctb"), person("Bert", "van Schaeybroeck", , "bertvs@meteo.be", role = "ctb")) -Description: Tools to exploit dynamical seasonal forecasts in order to provide information relevant to stakeholders at the seasonal timescale. The package contains process-based methods for forecast calibration, bias correction, statistical and stochastic downscaling, optimal forecast combination and multivariate verification, as well as basic and advanced tools to obtain tailored products. +Description: Tools to exploit dynamical seasonal forecasts in order to provide + information relevant to stakeholders at the seasonal timescale. The package + contains process-based methods for forecast calibration, bias correction, + statistical and stochastic downscaling, optimal forecast combination and + multivariate verification, as well as basic and advanced tools to obtain + tailored products. Depends: R (>= 3.2.0), maps Imports: s2dverification, - abind, - stats, - utils, - grDevices, - graphics, - JuliaCall, + rainfarmr, + multiApply, ncdf4, - ggplot2, plyr, + abind, data.table, reshape2, - multiApply + ggplot2, + graphics, + grDevices, + stats, + utils Suggests: zeallot, testthat, @@ -40,4 +46,4 @@ VignetteBuilder: R.rsp License: Apache License 2.0 | file LICENSE Encoding: UTF-8 LazyData: true -RoxygenNote: 6.1.1 +RoxygenNote: 5.0.0 diff --git a/NAMESPACE b/NAMESPACE index 872506f75a52612871a1f7a21a1cd06876fec13f..dae378533109f6ef0c5aee56cfe2b94deccfc8b8 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -7,16 +7,16 @@ export(CST_Load) export(CST_MultiMetric) export(CST_MultivarRMSE) export(CST_RFSlope) +export(CST_RFWeights) export(CST_RainFARM) export(PlotForecastPDF) export(PlotMostLikelyQuantileMap) export(RFSlope) -export(RFWeights) export(RainFARM) -import(JuliaCall) import(ggplot2) import(multiApply) import(ncdf4) +import(rainfarmr) import(s2dverification) import(stats) importFrom(data.table,CJ) @@ -45,6 +45,4 @@ importFrom(maps,map) importFrom(plyr,.) importFrom(plyr,dlply) importFrom(reshape2,melt) -importFrom(s2dverification,Load) importFrom(utils,glob2rx) -importFrom(utils,installed.packages) diff --git a/R/CST_Load.R b/R/CST_Load.R index 5e1dc3aaa5f20321e0ce148fd2f41c527438ca2a..66535583735c120572a770f6fc7cf9bbc50f017f 100644 --- a/R/CST_Load.R +++ b/R/CST_Load.R @@ -9,7 +9,7 @@ #' @param ... Parameters that are automatically forwarded to the `s2dverification::Load` function. See details in `?s2dverification::Load`. #' @return A list with one or two S3 objects, named 'exp' and 'obs', of the class 's2dv_cube', containing experimental and date-corresponding observational data, respectively. These 's2dv_cube's can be ingested by other functions in CSTools. If the parameter `exp` in the call to `CST_Load` is set to `NULL`, then only the 'obs' component is returned, and viceversa. #' @author Nicolau Manubens, \email{nicolau.manubens@bsc.es} -#' @importFrom s2dverification Load +#' @import s2dverification #' @importFrom utils glob2rx #' @export #' @examples @@ -36,7 +36,7 @@ #' obs <- CSTools::lonlat_data$obs #' } CST_Load <- function(...) { - exp <- s2dverification::Load(...) + exp <- Load(...) if (is.null(exp) || (is.null(exp$mod) && is.null(exp$obs))) { stop("The s2dverification::Load call did not return any data.") diff --git a/R/CST_MultivarRMSE.R b/R/CST_MultivarRMSE.R index 88967a134aca76502e3eb3bf1a32d18b9bf48217..c5ab53b5291a690d50a139b317db39b77b3ccad7 100644 --- a/R/CST_MultivarRMSE.R +++ b/R/CST_MultivarRMSE.R @@ -111,8 +111,7 @@ CST_MultivarRMSE <- function(exp, obs, weight = NULL) { AvgExp <- MeanListDim(exp[[j]]$data, narm = TRUE, c(2, 4)) AvgObs <- MeanListDim(obs[[j]]$data, narm = TRUE, c(2, 4)) # multivariate RMSE (weighted) - rmse <- s2dverification::RMS(AvgExp, AvgObs, posloop = 1, posRMS = 2, - conf = FALSE) + rmse <- RMS(AvgExp, AvgObs, posloop = 1, posRMS = 2, conf = FALSE) stdev <- sd(AvgObs) mvrmse <- mvrmse + (rmse / stdev * as.numeric(weight[j])) sumweights <- sumweights + as.numeric(weight[j]) diff --git a/R/CST_RFSlope.R b/R/CST_RFSlope.R index 6d7ae8af6a6b001a1a1cd68f0880b6fbc772aadc..6c996fc86df1659bb5cec4a75857333eb698e993 100644 --- a/R/CST_RFSlope.R +++ b/R/CST_RFSlope.R @@ -5,8 +5,6 @@ #' #' @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 (of the class 's2dv_cube') as 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 An object of the class 's2dv_cube', containing the spatial precipitation fields to downscale. #' The data object is expected to have an element named \code{$data} with at least two @@ -21,7 +19,7 @@ #' (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}. -#' importFrom utils installed.packages +#' @import rainfarmr #' @export #' @examples #' #Example using CST_RFSlope for a CSTools object @@ -59,8 +57,6 @@ CST_RFSlope <- function(data, kmin = 1, time_dim = NULL) { #' #' @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 @@ -80,6 +76,7 @@ CST_RFSlope <- function(data, kmin = 1, time_dim = NULL) { #' (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}. +#' @import multiApply #' @export #' @examples #' # Example for the 'reduced' RFSlope function @@ -139,20 +136,10 @@ RFSlope <- function(data, kmin = 1, time_dim = NULL, 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))) + data <- .aperm2(data, c(which(!imask), which(imask))) cdim <- dim(data) ind <- 1:length(which(!imask)) # compact (multiply) time_dim dimensions @@ -174,9 +161,8 @@ RFSlope <- function(data, kmin = 1, time_dim = NULL, .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") + fxp <- fft2d(pr) + sx <- fitslope(fxp, kmin = kmin) return(sx) } @@ -185,7 +171,7 @@ RFSlope <- function(data, kmin = 1, time_dim = NULL, # 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 ) + idim <- which(names(dim(field)) %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))) diff --git a/R/CST_RFWeights.R b/R/CST_RFWeights.R new file mode 100644 index 0000000000000000000000000000000000000000..03c828394bc4e5022848ac048ac94344997a5726 --- /dev/null +++ b/R/CST_RFWeights.R @@ -0,0 +1,84 @@ +#' 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. +#' @references 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 . +#' @param climfile 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 of longitudes. +#' @param lat Vector of latitudes. +#' The number of longitudes and latitudes is expected to be even and the same. If not +#' the function will perform a subsetting to ensure this condition. +#' @param varname Name of the variable to be read from \code{reffile}. +#' @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 rainfarmr +#' @examples +#' # Create weights to be used with the CST_RainFARM() or RainFARM() functions +#' # using an external fine-scale climatology file. +#' +#' \dontrun{ +#' # Specify lon and lat of the input +#' lon <- seq(10,13.5,0.5) +#' lat <- seq(40,43.5,0.5) +#' nf <- 8 +#' ww <- CST_RFWeights("./worldclim.nc", nf, lon, lat, fsmooth = TRUE) +#' } +#' @export + +CST_RFWeights <- function(climfile, nf, lon, lat, varname = "", fsmooth=TRUE) { + + # Ensure input grid is square and with even dimensions + if ((length(lat) != length(lon)) | (length(lon) %% 2 == 1)) { + print("Warning: input data are expected to be on a square grid") + print(" with an even number of pixels per side") + nmin <- min(length(lon), length(lat)) + nmin <- floor(nmin / 2) * 2 + lon <- lon[1:nmin] + lat <- lat[1:nmin] + print("The input data have been cut to the range") + print(paste0("lon: [", lon[1], ", ", lon[length(lon)], "] ", + " lat: [", lat[1], ", ", lat[length(lat)], "]")) + } + + ncin <- nc_open(climfile) + latin <- ncvar_get(ncin, grep("lat", attributes(ncin$dim)$names, value = TRUE)) + lonin <- ncvar_get(ncin, grep("lon", attributes(ncin$dim)$names, value = TRUE)) + if ( varname == "") { + varname <- grep("bnds", attributes(ncin$var)$names, + invert = TRUE, value = TRUE)[1] + } + zclim <- ncvar_get(ncin, varname) + +# Check if lon and lat need to be reversed + if ( lat[1] > lat[2] ) { + lat <- rev(lat) + frev = TRUE + } else { + frev = FALSE + } + if ( latin[1] > latin[2] ) { + latin <- rev(latin) + zclim = zclim[,seq(dim(zclim)[2],1)] + } + ww <- rfweights(zclim, lonin, latin, lon, lat, nf, fsmooth = fsmooth) + if ( frev ) { + ww = ww[,seq(dim(ww)[2],1)] + } + attributes(dim(ww))$names <- c("lon","lat") + return(ww) +} diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R index 1544a1e2ef314adae6ee941abaa24913df16c892..36800e55afb455fb0fa03c022219eccde1406f1c 100644 --- a/R/CST_RainFARM.R +++ b/R/CST_RainFARM.R @@ -8,12 +8,9 @@ #' 's2dv_cube' as provided by `CST_Load`) as 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 , +#' @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 An object of the class 's2dv_cube' as returned by `CST_Load`, #' containing the spatial precipitation fields to downscale. #' The data object is expected to have an element named \code{$data} with at least two @@ -23,7 +20,7 @@ #' The number of longitudes and latitudes in the input data is expected to be even and the same. If not #' the function will perform a subsetting to ensure this condition. #' @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 \code{CST_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). @@ -39,7 +36,7 @@ #' 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 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: #' @@ -57,7 +54,7 @@ #' \code{drop_realization_dim=TRUE} is specified). #' The ordering of the remaining dimensions in the \code{$data} element of the input object is maintained. #' @import multiApply -#' @importFrom utils installed.packages +#' @import rainfarmr #' @export #' @examples #' \donttest{ @@ -110,8 +107,6 @@ CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1, #' 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) @@ -122,7 +117,7 @@ CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1, #' @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 \code{CST_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). @@ -140,7 +135,7 @@ CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1, #' 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 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: #' @@ -156,7 +151,6 @@ CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1, #' 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 @@ -174,8 +168,10 @@ CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1, #' 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 -#' ww <- RFWeights("./worldclim.nc", nf, lon = lon_mat, lat = lat_mat, -#' weightsfn = "", fsmooth = TRUE) +#' \dontrun{ +#' ww <- CST_RFWeights("./worldclim.nc", nf, lon = lon_mat, lat = lat_mat, +#' fsmooth = TRUE) +#' } #' # downscale using weights (ww=1. means do not use weights) #' res <- RainFARM(pr, lon_mat, lat_mat, nf, #' fsmooth = TRUE, fglob = FALSE, @@ -242,25 +238,15 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1, 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]] + r <- lon_lat_fine(lon, lat, nf) + lon_f <- r$lon + lat_f <- r$lat # 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))) + data <- .aperm2(data, c(which(!imask), which(imask))) cdim <- dim(data) ind <- 1:length(which(!imask)) # compact (multiply) time_dim dimensions @@ -306,7 +292,7 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1, #' 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). +#' the \code{CST_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}. @@ -314,16 +300,15 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 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). +#' @param verbose Logical for verbose output (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") + fxp <- fft2d(pr) + sx <- fitslope(fxp, kmin = kmin) } else { sx <- slope } @@ -333,9 +318,8 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1, 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") + r[, , , i] <- rainfarm(pr, sx, nf, weights, fglob = fglob, + fsmooth = fsmooth, verbose = verbose) } return(r) } diff --git a/R/RFWeights.R b/R/RFWeights.R deleted file mode 100644 index 0459e285029489d210a4800ae491f27b082342a2..0000000000000000000000000000000000000000 --- a/R/RFWeights.R +++ /dev/null @@ -1,91 +0,0 @@ -#' 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 of longitudes (can be omitted if \code{reffile} is used). -#' @param lat Vector of latitudes (can be omitted if \code{reffile} is used). -#' The number of longitudes and latitudes is expected to be even and the same. If not -#' the function will perform a subsetting to ensure this condition. -#' @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) -#' 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 == "") { - # Ensure input grid is square and with even dimensions - if ((length(lat) != length(lon)) | (length(lon) %% 2 == 1)) { - print("Warning: input data are expected to be on a square grid") - print(" with an even number of pixels per side") - nmin <- min(length(lon), length(lat)) - nmin <- floor(nmin / 2) * 2 - lon <- lon[1:nmin] - lat <- lat[1:nmin] - print("The input data have been cut to the range") - print(paste0("lon: [", lon[1], ", ", lon[length(lon)], "] ", - " lat: [", lat[1], ", ", lat[length(lat)], "]")) - } - 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/R/sample_data.R b/R/sample_data.R index 3ceb7eedee75cd3beb1a7ae692b52866ea38c683..87ce3d0e4f7376d0e7cbc59f37515607724c4bbe 100644 --- a/R/sample_data.R +++ b/R/sample_data.R @@ -3,20 +3,20 @@ #' This sample data set contains gridded seasonal forecast and corresponding observational data from the Copernicus Climate Change ECMWF-System 5 forecast system, and from the Copernicus Climate Change ERA-5 reconstruction. Specifically, for the 'tas' (2-meter temperature) variable, for the 15 first forecast ensemble members, monthly averaged, for the 3 first forecast time steps (lead months 1 to 4) of the November start dates of 2000 to 2005, for the Mediterranean region (27N-48N, 12W-40E). The data was generated on (or interpolated onto, for the reconstruction) a rectangular regular grid of size 360 by 181. #' #' It is recommended to use the data set as follows: -#' ``` +#'\code{ #' require(zeallot) #' c(exp, obs) %<-% CSTools::lonlat_data -#' ``` +#'} #' #' The `CST_Load` call used to generate the data set in the infrastructure of the Earth Sciences Department of the Barcelona Supercomputing Center is shown next. Note that `CST_Load` internally calls `s2dverification::Load`, which would require a configuration file (not provided here) expressing the distribution of the 'system5c3s' and 'era5' NetCDF files in the file system. -#' ``` +#'\code{ #' library(CSTools) #' require(zeallot) #' #' startDates <- c('20001101', '20011101', '20021101', #' '20031101', '20041101', '20051101') #' -#' c(exp, obs) %<-% +#' lonlat_data <- #' CST_Load( #' var = 'tas', #' exp = 'system5c3s', @@ -29,7 +29,7 @@ #' output = 'lonlat', #' nprocs = 1 #' ) -#' ``` +#'} #' #' @name lonlat_data #' @docType data @@ -37,25 +37,52 @@ #' @keywords data NULL +#' Sample Of Experimental Precipitation Data In Function Of Longitudes And Latitudes +#' +#' This sample data set contains a small cutout of gridded seasonal precipitation forecast data from the Copernicus Climate Change ECMWF-System 5 forecast system, to be used to demonstrate downscaling. Specifically, for the 'pr' (precipitation) variable, for the first 6 forecast ensemble members, daily values, for all 31 days in March following the forecast starting dates in November of years 2010 to 2012, for a small 4x4 pixel cutout in a region in the North-Western Italian Alps (44N-47N, 6E-9E). The data resolution is 1 degree. +#' +#' The `CST_Load` call used to generate the data set in the infrastructure of the Marconi machine at CINECA is shown next, working on files which were extracted from forecast data available in the MEDSCOPE internal archive. +#' +#'\code{ +#' library(CSTools) +#' +#' infile <- list(path = '../medscope/nwalps/data/$VAR_NAME$_$START_DATE$_nwalps.nc') +#' lonlat_prec <- CST_Load('prlr', exp = list(infile), obs = NULL, +#' sdates = c('20101101', '20111101', '20121101'), +#' leadtimemin = 121, leadtimemax = 151, +#' latmin = 44, latmax = 47, +#' lonmin = 5, lonmax = 9, +#' nmember = 25, +#' storefreq = "daily", sampleperiod = 1, +#' output = "lonlat" +#' )$exp +#'} +#' +#' @name lonlat_prec +#' @docType data +#' @author Jost von Hardenberg \email{j.vonhardenberg@isac.cnr.it} +#' @keywords data +NULL + #' Sample Of Experimental And Observational Climate Data Averaged Over A Region #' #' This sample data set contains area-averaged seasonal forecast and corresponding observational data from the Copernicus Climate Change ECMWF-System 5 forecast system, and from the Copernicus Climate Change ERA-5 reconstruction. Specifically, for the 'tas' (2-meter temperature) variable, for the 15 first forecast ensemble members, monthly averaged, for the 3 first forecast time steps (lead months 1 to 4) of the November start dates of 2000 to 2005, for the Mediterranean region (27N-48N, 12W-40E). #' #' It is recommended to use the data set as follows: -#' ``` +#'\code{ #' require(zeallot) #' c(exp, obs) %<-% CSTools::areave_data -#' ``` +#'} #' #' The `CST_Load` call used to generate the data set in the infrastructure of the Earth Sciences Department of the Barcelona Supercomputing Center is shown next. Note that `CST_Load` internally calls `s2dverification::Load`, which would require a configuration file (not provided here) expressing the distribution of the 'system5c3s' and 'era5' NetCDF files in the file system. -#' ``` +#'\code{ #' library(CSTools) #' require(zeallot) #' #' startDates <- c('20001101', '20011101', '20021101', #' '20031101', '20041101', '20051101') #' -#' c(exp, obs) %<-% +#' areave_data <- #' CST_Load( #' var = 'tas', #' exp = 'system5c3s', @@ -68,7 +95,7 @@ NULL #' output = 'areave', #' nprocs = 1 #' ) -#' ``` +#'} #' #' @name areave_data #' @docType data diff --git a/R/zzz.R b/R/zzz.R new file mode 100644 index 0000000000000000000000000000000000000000..662624d9cb94cff6436d71016703c558c974c579 --- /dev/null +++ b/R/zzz.R @@ -0,0 +1,18 @@ +# Function to permute arrays of non-atomic elements (e.g. POSIXct) +.aperm2 <- function(x, new_order) { + old_dims <- dim(x) + attr_bk <- attributes(x) + if ('dim' %in% names(attr_bk)) { + attr_bk[['dim']] <- NULL + } + if (is.numeric(x)) { + x <- aperm(x, new_order) + } else { + y <- array(1:length(x), dim = dim(x)) + y <- aperm(y, new_order) + x <- x[as.vector(y)] + } + dim(x) <- old_dims[new_order] + attributes(x) <- c(attributes(x), attr_bk) + x +} diff --git a/data/lonlat_prec.RData b/data/lonlat_prec.RData new file mode 100644 index 0000000000000000000000000000000000000000..6725478b7a5c9c77613ad4c97e17d698f29e2402 Binary files /dev/null and b/data/lonlat_prec.RData differ diff --git a/man/CST_Anomaly.Rd b/man/CST_Anomaly.Rd index b214a31a299c6736e7704bec16378ba86ef7b0d4..e1c31f0ce5d3f09935f16adfa3c1a1fb08e620ca 100644 --- a/man/CST_Anomaly.Rd +++ b/man/CST_Anomaly.Rd @@ -53,12 +53,13 @@ str(anom3) anom4 <- CST_Anomaly(exp = exp, obs = obs, cross = FALSE, memb = FALSE) str(anom4) -} -\seealso{ -\code{\link[s2dverification]{Ano_CrossValid}}, \code{\link[s2dverification]{Clim}} and \code{\link{CST_Load}} } \author{ Perez-Zanon Nuria, \email{nuria.perez@bsc.es} Pena Jesus, \email{jesus.pena@bsc.es} } +\seealso{ +\code{\link[s2dverification]{Ano_CrossValid}}, \code{\link[s2dverification]{Clim}} and \code{\link{CST_Load}} +} + diff --git a/man/CST_BiasCorrection.Rd b/man/CST_BiasCorrection.Rd index a1b415fb525a9f2e3c72171e631c2633e0aefcd6..e8a82af0b9e67719a600e1e6b04e9fec2c0d96b6 100644 --- a/man/CST_BiasCorrection.Rd +++ b/man/CST_BiasCorrection.Rd @@ -35,9 +35,10 @@ attr(obs, 'class') <- 's2dv_cube' a <- CST_BiasCorrection(exp = exp, obs = obs) str(a) } -\references{ -Torralba, V., F.J. Doblas-Reyes, D. MacLeod, I. Christel and M. Davis (2017). Seasonal climate prediction: a new source of information for the management of wind energy resources. Journal of Applied Meteorology and Climatology, 56, 1231-1247, doi:10.1175/JAMC-D-16-0204.1. (CLIM4ENERGY, EUPORIAS, NEWA, RESILIENCE, SPECS) -} \author{ Verónica Torralba, \email{veronica.torralba@bsc.es} } +\references{ +Torralba, V., F.J. Doblas-Reyes, D. MacLeod, I. Christel and M. Davis (2017). Seasonal climate prediction: a new source of information for the management of wind energy resources. Journal of Applied Meteorology and Climatology, 56, 1231-1247, doi:10.1175/JAMC-D-16-0204.1. (CLIM4ENERGY, EUPORIAS, NEWA, RESILIENCE, SPECS) +} + diff --git a/man/CST_Calibration.Rd b/man/CST_Calibration.Rd index e3e402a5a0c4eeea83c2c6e3799da290dd4539a0..2e22da028d31be6c091331b5038baa15a83b1158 100644 --- a/man/CST_Calibration.Rd +++ b/man/CST_Calibration.Rd @@ -27,12 +27,13 @@ exp_calibrated <- CST_Calibration(exp = exp, obs = obs) str(exp_calibrated) } } +\author{ +Verónica Torralba, \email{veronica.torralba@bsc.es} +} \references{ Doblas-Reyes F.J, Hagedorn R, Palmer T.N. The rationale behind the success of multi-model ensembles in seasonal forecasting—II calibration and combination. Tellus A. 2005;57:234–252. doi:10.1111/j.1600-0870.2005.00104.x } \seealso{ \code{\link{CST_Load}} } -\author{ -Verónica Torralba, \email{veronica.torralba@bsc.es} -} + diff --git a/man/CST_Load.Rd b/man/CST_Load.Rd index bf03ba42d98810643a98bfca27f4563f6361162e..1fee022c23b928c072798a4136d927cbf6c2a287 100644 --- a/man/CST_Load.Rd +++ b/man/CST_Load.Rd @@ -47,3 +47,4 @@ obs <- CSTools::lonlat_data$obs \author{ Nicolau Manubens, \email{nicolau.manubens@bsc.es} } + diff --git a/man/CST_MultiMetric.Rd b/man/CST_MultiMetric.Rd index 99fdb400f1a82dd7267cd284a2b6597dfb314797..b6daf29ff29d9ed1cda38bd57101eb76a3a56cfa 100644 --- a/man/CST_MultiMetric.Rd +++ b/man/CST_MultiMetric.Rd @@ -37,14 +37,15 @@ c(ano_exp, ano_obs) \%<-\% CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb a <- CST_MultiMetric(exp = ano_exp, obs = ano_obs) str(a) } +\author{ +Mishra Niti, \email{niti.mishra@bsc.es} + +Perez-Zanon Nuria, \email{nuria.perez@bsc.es} +} \references{ Mishra, N., Prodhomme, C., & Guemas, V. (n.d.). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29–31.\url{http://link.springer.com/10.1007/s00382-018-4404-z} } \seealso{ \code{\link[s2dverification]{Corr}}, \code{\link[s2dverification]{RMS}}, \code{\link[s2dverification]{RMSSS}} and \code{\link{CST_Load}} } -\author{ -Mishra Niti, \email{niti.mishra@bsc.es} -Perez-Zanon Nuria, \email{nuria.perez@bsc.es} -} diff --git a/man/CST_MultivarRMSE.Rd b/man/CST_MultivarRMSE.Rd index 24af608c8360985de73763434a16f1a11fd35ffe..685eaf7712ccabbcef7d6cdea64a2f6f23556a31 100644 --- a/man/CST_MultivarRMSE.Rd +++ b/man/CST_MultivarRMSE.Rd @@ -56,9 +56,10 @@ weight <- c(1, 2) a <- CST_MultivarRMSE(exp = ano_exp, obs = ano_obs, weight = weight) str(a) } -\seealso{ -\code{\link[s2dverification]{RMS}} and \code{\link{CST_Load}} -} \author{ Deborah Verfaillie, \email{deborah.verfaillie@bsc.es} } +\seealso{ +\code{\link[s2dverification]{RMS}} and \code{\link{CST_Load}} +} + diff --git a/man/CST_RFSlope.Rd b/man/CST_RFSlope.Rd index 489e7a97ed29f57e308020e71cc0e330dba09761..d87192095e9333cd3c40031a60f1cc27ce071c2c 100644 --- a/man/CST_RFSlope.Rd +++ b/man/CST_RFSlope.Rd @@ -24,13 +24,10 @@ 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}. -importFrom utils installed.packages } \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 (of the class 's2dv_cube') as 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 @@ -56,3 +53,4 @@ slopes \author{ Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} } + diff --git a/man/RFWeights.Rd b/man/CST_RFWeights.Rd similarity index 62% rename from man/RFWeights.Rd rename to man/CST_RFWeights.Rd index 073cda877ce3926ef14d5ec03a0dfb396ec97186..d9bb48694fe3574665cb417ad52d5ed031512ef4 100644 --- a/man/RFWeights.Rd +++ b/man/CST_RFWeights.Rd @@ -1,14 +1,13 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RFWeights.R -\name{RFWeights} -\alias{RFWeights} +% Please edit documentation in R/CST_RFWeights.R +\name{CST_RFWeights} +\alias{CST_RFWeights} \title{Compute climatological weights for RainFARM stochastic precipitation downscaling} \usage{ -RFWeights(orofile, nf, lon = NULL, lat = NULL, reffile = "", - weightsfn = "", varname = "", fsmooth = TRUE) +CST_RFWeights(climfile, nf, lon, lat, varname = "", fsmooth = TRUE) } \arguments{ -\item{orofile}{Filename of a fine-scale precipitation climatology. +\item{climfile}{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. @@ -20,17 +19,12 @@ websites. The latter data will need to be converted to NetCDF format before bein \item{nf}{Refinement factor for downscaling (the output resolution is increased by this factor).} -\item{lon}{Vector of longitudes (can be omitted if \code{reffile} is used).} +\item{lon}{Vector of longitudes.} -\item{lat}{Vector of latitudes (can be omitted if \code{reffile} is used). +\item{lat}{Vector of latitudes. The number of longitudes and latitudes is expected to be even and the same. If not the function will perform a subsetting to ensure this condition.} -\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.} @@ -40,31 +34,27 @@ 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{ +\dontrun{ # Specify lon and lat of the input -lon_mat <- seq(10,13.5,0.5) -lat_mat <- seq(40,43.5,0.5) - +lon <- seq(10,13.5,0.5) +lat <- 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") +ww <- CST_RFWeights("./worldclim.nc", nf, lon, lat, fsmooth = TRUE) } } \author{ Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} } +\references{ +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 . +} + diff --git a/man/CST_RainFARM.Rd b/man/CST_RainFARM.Rd index 5343af7aaf8599ecd997cdb300576f2ee1fcf100..e6731b2b0fb0468036123858a10301c0592a3a26 100644 --- a/man/CST_RainFARM.Rd +++ b/man/CST_RainFARM.Rd @@ -21,7 +21,7 @@ the function will perform a subsetting to ensure this condition.} \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 \code{CST_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.} @@ -43,7 +43,7 @@ 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{verbose}{Logical for verbose output (default: FALSE).} \item{drop_realization_dim}{Logical to remove the "realization" stochastic ensemble dimension (default: FALSE) with the following behaviour if set to TRUE: @@ -69,11 +69,6 @@ downscaling method and accepts a CSTools object (an object of the class 's2dv_cube' as provided by `CST_Load`) as 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{ \donttest{ @@ -102,3 +97,9 @@ dim(res$data) \author{ Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} } +\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. +} + diff --git a/man/PlotForecastPDF.Rd b/man/PlotForecastPDF.Rd index 8b8bae0f2feaa99cbe3e855bc7c81a9d526df6cc..562827ee44897010b46a41065661f575295e58ba 100644 --- a/man/PlotForecastPDF.Rd +++ b/man/PlotForecastPDF.Rd @@ -4,11 +4,10 @@ \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.} @@ -50,3 +49,4 @@ PlotForecastPDF(fcsts2, c(-0.66, 0.66), extreme.limits = c(-1.2, 1.2), \author{ Llorenç Lledó \email{llledo@bsc.es} } + diff --git a/man/PlotMostLikelyQuantileMap.Rd b/man/PlotMostLikelyQuantileMap.Rd index 0d984ede81188274f1957f6badf3c8fdf2add781..6c92850ead6ac72be3985fc233a273a7969e4060 100644 --- a/man/PlotMostLikelyQuantileMap.Rd +++ b/man/PlotMostLikelyQuantileMap.Rd @@ -109,10 +109,11 @@ PlotMostLikelyQuantileMap(bins, lons, lats, mask = 1 - (w1 + w2 / max(c(w1, w2))), brks = 20, width = 10, height = 8) -} -\seealso{ -\code{PlotCombinedMap} and \code{PlotEquiMap} } \author{ Veronica Torralba, \email{veronica.torralba@bsc.es}, Nicolau Manubens, \email{nicolau.manubens@bsc.es} } +\seealso{ +\code{PlotCombinedMap} and \code{PlotEquiMap} +} + diff --git a/man/RFSlope.Rd b/man/RFSlope.Rd index 628838b7b84aaa7de866c54a5c5671f2ff4af5fe..c4864e15ddd9cc11a02bf8fc60f95f0e7a4043c3 100644 --- a/man/RFSlope.Rd +++ b/man/RFSlope.Rd @@ -36,8 +36,6 @@ minus the dimensions specified by \code{lon_dim}, \code{lat_dim} and \code{time_ \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 @@ -65,3 +63,4 @@ slopes \author{ Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} } + diff --git a/man/RainFARM.Rd b/man/RainFARM.Rd index 34a4d418fcfb50c9bec63fc427c46177feb8b212..c0bd5f7dc5a2fa97eb3e6f294c6c1a2576e32845 100644 --- a/man/RainFARM.Rd +++ b/man/RainFARM.Rd @@ -4,10 +4,9 @@ \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) +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. @@ -25,7 +24,7 @@ the function will perform a subsetting to ensure this condition.} \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 \code{CST_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.} @@ -60,7 +59,7 @@ is searched and the first one with more than one element is chosen.} 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).} +\item{verbose}{logical for verbose output (default: FALSE).} } \value{ RainFARM() returns a list containing the fine-scale longitudes, latitudes @@ -78,8 +77,6 @@ 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{ \donttest{ @@ -96,8 +93,10 @@ lat_mat <- seq(40, 43.5, 0.5) 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 -ww <- RFWeights("./worldclim.nc", nf, lon = lon_mat, lat = lat_mat, - weightsfn = "", fsmooth = TRUE) +\dontrun{ +ww <- CST_RFWeights("./worldclim.nc", nf, lon = lon_mat, lat = lat_mat, + fsmooth = TRUE) +} # downscale using weights (ww=1. means do not use weights) res <- RainFARM(pr, lon_mat, lat_mat, nf, fsmooth = TRUE, fglob = FALSE, @@ -116,3 +115,4 @@ dim(res$data) \author{ Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} } + diff --git a/man/areave_data.Rd b/man/areave_data.Rd index d42663944e89e94abe084df648badc329e2778de..fd0f5abb41af01e0c95818ca8a1ed0f308974460 100644 --- a/man/areave_data.Rd +++ b/man/areave_data.Rd @@ -9,20 +9,20 @@ This sample data set contains area-averaged seasonal forecast and corresponding } \details{ It is recommended to use the data set as follows: -``` +\code{ require(zeallot) c(exp, obs) %<-% CSTools::areave_data -``` +} The `CST_Load` call used to generate the data set in the infrastructure of the Earth Sciences Department of the Barcelona Supercomputing Center is shown next. Note that `CST_Load` internally calls `s2dverification::Load`, which would require a configuration file (not provided here) expressing the distribution of the 'system5c3s' and 'era5' NetCDF files in the file system. -``` +\code{ library(CSTools) require(zeallot) startDates <- c('20001101', '20011101', '20021101', '20031101', '20041101', '20051101') -c(exp, obs) %<-% +areave_data <- CST_Load( var = 'tas', exp = 'system5c3s', @@ -35,9 +35,10 @@ c(exp, obs) %<-% output = 'areave', nprocs = 1 ) -``` +} } \author{ Nicolau Manubens \email{nicolau.manubens@bsc.es} } \keyword{data} + diff --git a/man/lonlat_data.Rd b/man/lonlat_data.Rd index a835ad89ef8890fc136265f7ee721481baaa81a6..8d0dbf1883d8aa6bae9eeba9afeef87d3b6fa6a9 100644 --- a/man/lonlat_data.Rd +++ b/man/lonlat_data.Rd @@ -9,20 +9,20 @@ This sample data set contains gridded seasonal forecast and corresponding observ } \details{ It is recommended to use the data set as follows: -``` +\code{ require(zeallot) c(exp, obs) %<-% CSTools::lonlat_data -``` +} The `CST_Load` call used to generate the data set in the infrastructure of the Earth Sciences Department of the Barcelona Supercomputing Center is shown next. Note that `CST_Load` internally calls `s2dverification::Load`, which would require a configuration file (not provided here) expressing the distribution of the 'system5c3s' and 'era5' NetCDF files in the file system. -``` +\code{ library(CSTools) require(zeallot) startDates <- c('20001101', '20011101', '20021101', '20031101', '20041101', '20051101') -c(exp, obs) %<-% +lonlat_data <- CST_Load( var = 'tas', exp = 'system5c3s', @@ -35,9 +35,10 @@ c(exp, obs) %<-% output = 'lonlat', nprocs = 1 ) -``` +} } \author{ Nicolau Manubens \email{nicolau.manubens@bsc.es} } \keyword{data} + diff --git a/man/lonlat_prec.Rd b/man/lonlat_prec.Rd new file mode 100644 index 0000000000000000000000000000000000000000..23678ead65dcde1a279c6adfe5a9142da95c5aa1 --- /dev/null +++ b/man/lonlat_prec.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/sample_data.R +\docType{data} +\name{lonlat_prec} +\alias{lonlat_prec} +\title{Sample Of Experimental Precipitation Data In Function Of Longitudes And Latitudes} +\description{ +This sample data set contains a small cutout of gridded seasonal precipitation forecast data from the Copernicus Climate Change ECMWF-System 5 forecast system, to be used to demonstrate downscaling. Specifically, for the 'pr' (precipitation) variable, for the first 6 forecast ensemble members, daily values, for all 31 days in March following the forecast starting dates in November of years 2010 to 2012, for a small 4x4 pixel cutout in a region in the North-Western Italian Alps (44N-47N, 6E-9E). The data resolution is 1 degree. +} +\details{ +The `CST_Load` call used to generate the data set in the infrastructure of the Marconi machine at CINECA is shown next, working on files which were extracted from forecast data available in the MEDSCOPE internal archive. + +\code{ +library(CSTools) + +infile <- list(path = '../medscope/nwalps/data/$VAR_NAME$_$START_DATE$_nwalps.nc') +lonlat_prec <- CST_Load('prlr', exp = list(infile), obs = NULL, + sdates = c('20101101', '20111101', '20121101'), + leadtimemin = 121, leadtimemax = 151, + latmin = 44, latmax = 47, + lonmin = 5, lonmax = 9, + nmember = 25, + storefreq = "daily", sampleperiod = 1, + output = "lonlat" + )$exp +} +} +\author{ +Jost von Hardenberg \email{j.vonhardenberg@isac.cnr.it} +} +\keyword{data} + diff --git a/tests/testthat/test-CST_BiasCorrection.R b/tests/testthat/test-CST_BiasCorrection.R index 39d5daa510a2c85e46daf49b7758bd60b5d055de..c5525c57012ec0b92c46e7bc9e2284831a26a66e 100644 --- a/tests/testthat/test-CST_BiasCorrection.R +++ b/tests/testthat/test-CST_BiasCorrection.R @@ -17,13 +17,14 @@ test_that("Sanity checks", { obs <- list(data = obs1, lat = lat, lon = lon) attr(exp, 'class') <- 's2dv_cube' attr(obs, 'class') <- 's2dv_cube' - - expect_equal(length(CST_BiasCorrection(exp = exp, obs = obs)), 3) - expect_equal(dim(CST_BiasCorrection(exp = exp, obs = obs)$data), + + bc <- CST_BiasCorrection(exp = exp, obs = obs) + expect_equal(length(bc), 3) + expect_equal(dim(bc$data), c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)) - expect_equal(CST_BiasCorrection(exp = exp, obs = obs)$lat, lat) - expect_equal(CST_BiasCorrection(exp = exp, obs = obs)$lon, lon) + expect_equal(bc$lat, lat) + expect_equal(bc$lon, lon) expect_error(CST_BiasCorrection(exp = exp, obs = exp), paste0("The length of the dimension 'member' in the component 'data' ", diff --git a/tests/testthat/test-CST_Calibration.R b/tests/testthat/test-CST_Calibration.R index 23cd4c21311c06e83b27c79c48a3856c0d8b5913..917b6ac5bd31c565709881b7ee7e1918c29d9d9f 100644 --- a/tests/testthat/test-CST_Calibration.R +++ b/tests/testthat/test-CST_Calibration.R @@ -2,34 +2,42 @@ context("Generic tests") test_that("Sanity checks", { expect_error( CST_Calibration(exp = 1), - c("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ", - "as output by CSTools::CST_Load.")) + "Parameter 'exp' and 'obs' must be of the class 's2dv_cube', " + ) expect_error( CST_Calibration(obs = 1), - c("argument \"exp\" is missing, with no default")) - + c("argument \"exp\" is missing, with no default") + ) library(zeallot) c(exp, obs) %<-% lonlat_data - expect_equal(length(CST_Calibration(exp = exp, obs = obs)), 9) - expect_equal(dim(CST_Calibration(exp = exp, obs = obs)$data), dim(exp$data)) - expect_equal(CST_Calibration(exp = exp, obs = obs)$lat, exp$lat) - expect_equal(CST_Calibration(exp = exp, obs = obs)$lat, obs$lat) - expect_equal(CST_Calibration(exp = exp, obs = obs)$lon, exp$lon) - expect_equal(CST_Calibration(exp = exp, obs = obs)$lon, obs$lon) - expect_error(CST_Calibration(exp = exp, obs = exp), - c("The length of the dimension 'member' in the component 'data' ", - "of the parameter 'obs' must be equal to 1.")) - + cal <- CST_Calibration(exp = exp, obs = obs) + expect_equal(length(cal), 9) + expect_equal(dim(cal$data), dim(exp$data)) + expect_equal(cal$lat, exp$lat) + expect_equal(cal$lat, obs$lat) + expect_equal(cal$lon, exp$lon) + expect_equal(cal$lon, obs$lon) + expect_error( + CST_Calibration(exp = exp, obs = exp), + "The length of the dimension 'member' in the component 'data' " + ) + exp2 <- exp exp2$data[1, 2, 1, 1, 1, 1] <- NA - expect_warning(CST_Calibration(exp = exp2, obs = obs), - "Parameter 'exp' contains NA values.") - + expect_warning( + CST_Calibration(exp = exp2, obs = obs), + "Parameter 'exp' contains NA values." + ) + obs2 <- obs obs2$data[1, 1, 2, 1, 1, 1] <- NA - expect_warning(CST_Calibration(exp = exp, obs = obs2), - "Parameter 'obs' contains NA values.") - - expect_warning(CST_Calibration(exp = exp2, obs = obs2), - "Parameter 'obs' contains NA values", "Parameter 'exp' contains NA values.") + expect_warning( + CST_Calibration(exp = exp, obs = obs2), + "Parameter 'obs' contains NA values." + ) + + expect_warning( + CST_Calibration(exp = exp2, obs = obs2), + "Parameter 'obs' contains NA values", "Parameter 'exp' contains NA values." + ) }) diff --git a/vignettes/RainFARM_vignette.md b/vignettes/RainFARM_vignette.md index f8c286436caa9452ff826818339f25a360b4c337..7b9817d19e5344bd474029d2e2f409e08d9c922c 100644 --- a/vignettes/RainFARM_vignette.md +++ b/vignettes/RainFARM_vignette.md @@ -3,12 +3,12 @@ title: "Rainfall Filtered Autoregressive Model (RainFARM) precipitation downscal author: "Jost von Hardenberg (ISAC-CNR)" date: "26/03/2019" output: - pdf_document: default html_document: default + pdf_document: default vignette: > - %\VignetteEngine{R.rsp::md} %\VignetteIndexEntry{RainFARM} %\usepackage[utf8]{inputenc} + %\VignetteEngine{R.rsp::md} --- -![The same field as in Fig. 1 downscaled using climatological weights from the WORLDCLIM dataset. (left panel). The center and right panel show the climatology of the downscaled fields over all ensemble members, realizations and forecast times without and with climatological weights respectively. The center and right panels use a different colorscale than the left panel.](Figures/RainFARM_fig2.png){width=640px} +![The same field as in Fig. 1 downscaled using climatological weights from the WORLDCLIM dataset. (left panel). The center and right panel show the climatology of the downscaled fields over all ensemble members, realizations and forecast times without and with climatological weights respectively, for the first starting date. The center and right panels use a different colorscale than the left panel.](Figures/RainFARM_fig2.png){width=640px} ### Determining the spectral slopes @@ -200,11 +196,6 @@ RainFARM creates an additional dimension in the `$data` array of the output data 3) if `nens > 1` and a "member" dimension does not exist: the "realization" dimension is renamed to "member". -## Technical issues - -For maximum efficiency, the `CST_RainFARM()`, `RainFARM()`, `RFWeights()` and `RFSlope()` functions are based on an external library written in the Julia language (https://github.com/jhardenberg/RainFARM.jl). The `JuliaCall` R package (https://CRAN.R-project.org/package=JuliaCall) is used to call the external library. -In particular the function uses the `julia_setup()` function from the JuliaCall library, which takes a long time (up to a few tens of seconds) only the very first time it is called. Further calls to these functions will present no delay. - ## Bibliography * Terzago, S., Palazzi, E., von Hardenberg, J. 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, doi:10.5194/nhess-18-2825-2018 (2018)