From 63a85d676d31d7c80076d92e43fb04da0306dc29 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <hardenberg@gmail.com>
Date: Sat, 27 Oct 2018 02:30:40 +0200
Subject: [PATCH 01/48] First working RainFARM function

---
 R/rainfarm.R | 46 ++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 46 insertions(+)
 create mode 100644 R/rainfarm.R

diff --git a/R/rainfarm.R b/R/rainfarm.R
new file mode 100644
index 00000000..b65c68e0
--- /dev/null
+++ b/R/rainfarm.R
@@ -0,0 +1,46 @@
+#'RainFARM Stochastic Precipitation Downscaling
+#'
+#'@author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
+#'@description This function implements the RainFARM stochastic precipitation downscaling method. Adapted for climate downscaling and including orographic correction. References: D'Onofrio et al. 2014, J of Hydrometeorology 15, 830-843; Rebora et. al 2006, JHM 7, 724
+#'Has still to be adapted for s2dverification
+#'@import JuliaCall
+#'@examples 
+#'# julia_setup takes a long time only the first time it is called on a machine
+#'
+#'library(JuliaCall)
+#'julia_setup()
+#'julia_library("RainFARM")
+#'source("rainfarm.R")
+#'
+#'nf=8; nens=3;
+#'# reads a test file with dimension 8x8 and 100 timesteps
+#'r=julia_call("read_netcdf2d","pr_ls.nc","pr", need_return = "R" );
+#'pr=r[[1]]; lon_mat=r[[2]]; lat_mat=r[[3]];
+#'ww <- array(1., dim=c(8*nf,8*nf));
+#'rd=rainfarm(pr,lon_mat,lat_mat,nf,fsmooth=TRUE,fglob=FALSE,weights=ww,nens=2,verbose=TRUE)
+#'@export
+
+rainfarm <- function(pr, lon_mat, lat_mat, nf, fglob=FALSE,fsmooth=TRUE,weights=1.,nens=1,verbose=FALSE,slope=0) {
+
+library(JuliaCall)
+julia_setup()
+julia_library("RainFARM")
+
+#r=julia_call("lon_lat_fine",lon_mat, lat_mat,nf,need_return = "R" );
+#lon_f=r[[1]]; lat_f=r[[2]];
+
+if(slope==0){
+r=julia_call("fft3d",pr,need_return = "R"  );
+fxp=r[[1]]; ftp=r[[2]];
+sx=julia_call("fitslopex",fxp,kmin=1,need_return = "R" );
+} else {
+sx=slope;
+}
+
+rd=list();
+for (j in 1:nens) {
+   rd=c(rd,list(julia_call("rainfarmn", pr, sx, nf, weights ,fglob=fglob,fsmooth=fsmooth,verbose=verbose, need_return = "R")));
+}
+
+return(rd)
+}
-- 
GitLab


From 0602efc43140be67410f45b5639f91edde75dd57 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <hardenberg@gmail.com>
Date: Sun, 28 Oct 2018 15:21:06 +0100
Subject: [PATCH 02/48] renamed rainfarmn

---
 R/rainfarm.R | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/R/rainfarm.R b/R/rainfarm.R
index b65c68e0..78b1487f 100644
--- a/R/rainfarm.R
+++ b/R/rainfarm.R
@@ -1,7 +1,9 @@
 #'RainFARM Stochastic Precipitation Downscaling
 #'
 #'@author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
-#'@description This function implements the RainFARM stochastic precipitation downscaling method. Adapted for climate downscaling and including orographic correction. References: D'Onofrio et al. 2014, J of Hydrometeorology 15, 830-843; Rebora et. al 2006, JHM 7, 724
+#'@description This function implements the RainFARM stochastic precipitation downscaling method. 
+#'Adapted for climate downscaling and including orographic correction. 
+#'References: D'Onofrio et al. 2014, J of Hydrometeorology 15, 830-843; Rebora et. al 2006, JHM 7, 724
 #'Has still to be adapted for s2dverification
 #'@import JuliaCall
 #'@examples 
@@ -26,8 +28,8 @@ library(JuliaCall)
 julia_setup()
 julia_library("RainFARM")
 
-#r=julia_call("lon_lat_fine",lon_mat, lat_mat,nf,need_return = "R" );
-#lon_f=r[[1]]; lat_f=r[[2]];
+r=julia_call("lon_lat_fine",lon_mat, lat_mat,nf,need_return = "R" );
+lon_f=r[[1]]; lat_f=r[[2]];
 
 if(slope==0){
 r=julia_call("fft3d",pr,need_return = "R"  );
@@ -39,8 +41,8 @@ sx=slope;
 
 rd=list();
 for (j in 1:nens) {
-   rd=c(rd,list(julia_call("rainfarmn", pr, sx, nf, weights ,fglob=fglob,fsmooth=fsmooth,verbose=verbose, need_return = "R")));
+   rd=c(rd,list(julia_call("rainfarm", pr, sx, nf, weights ,fglob=fglob,fsmooth=fsmooth,verbose=verbose, need_return = "R")));
 }
 
-return(rd)
+return(c(list(lon_f),list(lat_f),list(rd)))
 }
-- 
GitLab


From 2430e6c9ebcb795e67df676f6b717173f70d57d3 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <jost@Air-di-Jost.lan>
Date: Thu, 13 Dec 2018 21:41:17 +0100
Subject: [PATCH 03/48] style fixes

---
 R/rainfarm.R | 88 ++++++++++++++++++++++++++++------------------------
 1 file changed, 47 insertions(+), 41 deletions(-)

diff --git a/R/rainfarm.R b/R/rainfarm.R
index 78b1487f..09fa1e6e 100644
--- a/R/rainfarm.R
+++ b/R/rainfarm.R
@@ -1,48 +1,54 @@
-#'RainFARM Stochastic Precipitation Downscaling
+#' RainFARM Stochastic Precipitation Downscaling
 #'
-#'@author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
-#'@description This function implements the RainFARM stochastic precipitation downscaling method. 
-#'Adapted for climate downscaling and including orographic correction. 
-#'References: D'Onofrio et al. 2014, J of Hydrometeorology 15, 830-843; Rebora et. al 2006, JHM 7, 724
-#'Has still to be adapted for s2dverification
-#'@import JuliaCall
-#'@examples 
-#'# julia_setup takes a long time only the first time it is called on a machine
-#'
-#'library(JuliaCall)
-#'julia_setup()
-#'julia_library("RainFARM")
-#'source("rainfarm.R")
-#'
-#'nf=8; nens=3;
-#'# reads a test file with dimension 8x8 and 100 timesteps
-#'r=julia_call("read_netcdf2d","pr_ls.nc","pr", need_return = "R" );
-#'pr=r[[1]]; lon_mat=r[[2]]; lat_mat=r[[3]];
-#'ww <- array(1., dim=c(8*nf,8*nf));
-#'rd=rainfarm(pr,lon_mat,lat_mat,nf,fsmooth=TRUE,fglob=FALSE,weights=ww,nens=2,verbose=TRUE)
-#'@export
+#' @author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
+#' @description This function implements the RainFARM stochastic precipitation downscaling method.
+#' Adapted for climate downscaling and including orographic correction.
+#' References: D'Onofrio et al. 2014, J of Hydrometeorology 15, 830-843; Rebora et. al 2006, JHM 7, 724
+#' Has still to be adapted for s2dverification
+#' @import JuliaCall
+#' @examples
+#' # julia_setup takes a long time only the first time it is called on a machine
+#' 
+#' library(JuliaCall)
+#' julia_setup()
+#' julia_library("RainFARM")
+#' source("rainfarm.R")
+#' 
+#' nf <- 8
+#' nens <- 3
+#' # reads a test file with dimension 8x8 and 100 timesteps
+#' r <- julia_call("read_netcdf2d", "pr_ls.nc", "pr", need_return = "R")
+#' pr <- r[[1]]
+#' lon_mat <- r[[2]]
+#' lat_mat <- r[[3]]
+#' ww <- array(1., dim = c(8 * nf, 8 * nf))
+#' rd <- rainfarm(pr, lon_mat, lat_mat, nf, fsmooth = TRUE, fglob = FALSE, weights = ww, nens = 2, verbose = TRUE)
+#' @export
 
-rainfarm <- function(pr, lon_mat, lat_mat, nf, fglob=FALSE,fsmooth=TRUE,weights=1.,nens=1,verbose=FALSE,slope=0) {
+rainfarm <- function(pr, lon_mat, lat_mat, nf, fglob = FALSE, fsmooth = TRUE,
+                     weights = 1., nens = 1, verbose = FALSE, slope = 0) {
+  library(JuliaCall)
+  julia_setup()
+  julia_library("RainFARM")
 
-library(JuliaCall)
-julia_setup()
-julia_library("RainFARM")
+  r <- julia_call("lon_lat_fine", lon_mat, lat_mat, nf, need_return = "R")
+  lon_f <- r[[1]]
+  lat_f <- r[[2]]
 
-r=julia_call("lon_lat_fine",lon_mat, lat_mat,nf,need_return = "R" );
-lon_f=r[[1]]; lat_f=r[[2]];
+  if (slope == 0) {
+    r <- julia_call("fft3d", pr, need_return = "R")
+    fxp <- r[[1]]
+    ftp <- r[[2]]
+    sx <- julia_call("fitslopex", fxp, kmin = 1, need_return = "R")
+  } else {
+    sx <- slope
+  }
 
-if(slope==0){
-r=julia_call("fft3d",pr,need_return = "R"  );
-fxp=r[[1]]; ftp=r[[2]];
-sx=julia_call("fitslopex",fxp,kmin=1,need_return = "R" );
-} else {
-sx=slope;
-}
-
-rd=list();
-for (j in 1:nens) {
-   rd=c(rd,list(julia_call("rainfarm", pr, sx, nf, weights ,fglob=fglob,fsmooth=fsmooth,verbose=verbose, need_return = "R")));
-}
+  rd <- list()
+  for (j in 1:nens) {
+    rd <- c(rd, list(julia_call("rainfarm", pr, sx, nf, weights, fglob = fglob,
+                     fsmooth = fsmooth, verbose = verbose, need_return = "R")))
+  }
 
-return(c(list(lon_f),list(lat_f),list(rd)))
+  return(c(list(lon_f), list(lat_f), list(rd)))
 }
-- 
GitLab


From f5dd3839f79a6e4a3b5bc39ab2180705b59d8a86 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Fri, 21 Dec 2018 08:43:40 +0100
Subject: [PATCH 04/48] Added rfweights function

---
 R/rainfarm.R | 90 ++++++++++++++++++++++++++++++++++++++++------------
 1 file changed, 69 insertions(+), 21 deletions(-)

diff --git a/R/rainfarm.R b/R/rainfarm.R
index 09fa1e6e..50b810a6 100644
--- a/R/rainfarm.R
+++ b/R/rainfarm.R
@@ -1,53 +1,101 @@
 #' RainFARM Stochastic Precipitation Downscaling
 #'
 #' @author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
+#'
 #' @description This function implements the RainFARM stochastic precipitation downscaling method.
 #' Adapted for climate downscaling and including orographic correction.
 #' References: D'Onofrio et al. 2014, J of Hydrometeorology 15, 830-843; Rebora et. al 2006, JHM 7, 724
-#' Has still to be adapted for s2dverification
+#' The functions use julia_setup from the JuliaCall library, 
+#' which takes a long time only the first time it is called on a machine
+#'
 #' @import JuliaCall
 #' @examples
-#' # julia_setup takes a long time only the first time it is called on a machine
 #' 
-#' library(JuliaCall)
-#' julia_setup()
-#' julia_library("RainFARM")
 #' source("rainfarm.R")
+#' source("rfweights.R")
 #' 
-#' nf <- 8
-#' nens <- 3
-#' # reads a test file with dimension 8x8 and 100 timesteps
-#' r <- julia_call("read_netcdf2d", "pr_ls.nc", "pr", need_return = "R")
-#' pr <- r[[1]]
-#' lon_mat <- r[[2]]
-#' lat_mat <- r[[3]]
+#' 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(8,8,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))
-#' rd <- rainfarm(pr, lon_mat, lat_mat, nf, fsmooth = TRUE, fglob = FALSE, weights = ww, nens = 2, verbose = TRUE)
+#' 
+#' # ... or create proper weights using an external fine-scale climatology file
+#' #     Specify a weightsfn filename if you wish to save the weights
+#' 
+#' rd <- rfweights("./worldclim.nc", nf, lon=lon_mat, lat=lat_mat, weightsfn="", fsmooth = TRUE)
+#' 
+#' # downscale using weights (ww=1. means do not use weights)
+#' r <- rainfarm(pr, lon_mat, lat_mat, nf, fsmooth = TRUE, fglob = FALSE, weights = ww, nens = 2, verbose = TRUE)
+#' rd <- r[[1]]
+#' lon_f <- r[[2]]
+#' lat_f <- r[[3]]
+#'
 #' @export
 
-rainfarm <- function(pr, lon_mat, lat_mat, nf, fglob = FALSE, fsmooth = TRUE,
-                     weights = 1., nens = 1, verbose = FALSE, slope = 0) {
+rfweights <- function(orofile, nf, lon=NULL, lat=NULL, reffile="",
+                      weightsfn="", varname="", fsmooth=TRUE) {
+
+  library(JuliaCall)
+  library(ncdf4)
+  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 <- "__reffile.nc"
+     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)
+
+  system("rm -f __reffile.nc gridrf_2_*.nc")
+  return(ww)
+}
+
+rainfarm <- function(pr, lon, lat, nf, fglob = FALSE,
+                     fsmooth = TRUE, weights = 1., nens = 1,
+                     verbose = FALSE, slope = 0, kmin=1) {
   library(JuliaCall)
-  julia_setup()
+  julia_setup(useRCall = FALSE)
   julia_library("RainFARM")
 
-  r <- julia_call("lon_lat_fine", lon_mat, lat_mat, nf, need_return = "R")
+  r <- julia_call("lon_lat_fine", lon, lat, as.integer(nf),
+                  need_return = "R")
   lon_f <- r[[1]]
   lat_f <- r[[2]]
 
   if (slope == 0) {
     r <- julia_call("fft3d", pr, need_return = "R")
     fxp <- r[[1]]
-    ftp <- r[[2]]
-    sx <- julia_call("fitslopex", fxp, kmin = 1, need_return = "R")
+    sx <- julia_call("fitslopex", fxp, kmin = as.integer(kmin),
+                     need_return = "R")
   } else {
     sx <- slope
   }
 
   rd <- list()
   for (j in 1:nens) {
-    rd <- c(rd, list(julia_call("rainfarm", pr, sx, nf, weights, fglob = fglob,
-                     fsmooth = fsmooth, verbose = verbose, need_return = "R")))
+    rd <- c(rd, list(julia_call("rainfarm", pr, sx, as.integer(nf), weights,
+                                fglob = fglob, fsmooth = fsmooth,
+                                verbose = verbose, need_return = "R")))
   }
 
   return(c(list(lon_f), list(lat_f), list(rd)))
-- 
GitLab


From 6305c29d7fd39a8a7c980f5996d742ff92303bd8 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Fri, 21 Dec 2018 09:00:44 +0100
Subject: [PATCH 05/48] updated imports and DESCRIPTION

---
 DESCRIPTION  | 6 +++++-
 R/rainfarm.R | 2 +-
 2 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index d525653c..be0678f5 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -10,6 +10,8 @@ Authors@R: c(
     person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut", "cre")),
     person("Deborah", "Verfaillie", , "deborah.verfaillie@bsc.es", role = "aut"),
     person("Nicolau", "Manubens", , "nicolau.manubens@bsc.es", role = "ctb"),
+    person("CNR-ISAC", role = c("aut", "cph")),
+    person("Jost", "von Hardenberg", "j.vonhardenberg@isac.cnr.it", role = c("aut", "cph")),
     person("Other authors", role = "aut"),
     person("Other contributors", role = "ctb"))
 Description: Improved global climate model calibration and regionalization
@@ -21,7 +23,9 @@ Depends:
 Imports:
     s2dverification,
     abind,
-    stats
+    stats,
+    JuliaCall,
+    ncdf4
 License: LGPL-3
 URL: https://earth.bsc.es/gitlab/external/medscope
 BugReports: https://earth.bsc.es/gitlab/external/medscope/issues
diff --git a/R/rainfarm.R b/R/rainfarm.R
index 50b810a6..d8693aee 100644
--- a/R/rainfarm.R
+++ b/R/rainfarm.R
@@ -8,7 +8,7 @@
 #' The functions use julia_setup from the JuliaCall library, 
 #' which takes a long time only the first time it is called on a machine
 #'
-#' @import JuliaCall
+#' @import JuliaCall, ncdf4
 #' @examples
 #' 
 #' source("rainfarm.R")
-- 
GitLab


From 4708330a1aea6488fd8922f6e21659734669961b Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Fri, 21 Dec 2018 09:04:21 +0100
Subject: [PATCH 06/48] typo in description

---
 R/rainfarm.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/rainfarm.R b/R/rainfarm.R
index d8693aee..3b9d2120 100644
--- a/R/rainfarm.R
+++ b/R/rainfarm.R
@@ -31,7 +31,7 @@
 #' # ... or create proper weights using an external fine-scale climatology file
 #' #     Specify a weightsfn filename if you wish to save the weights
 #' 
-#' rd <- rfweights("./worldclim.nc", nf, lon=lon_mat, lat=lat_mat, weightsfn="", fsmooth = TRUE)
+#' ww <- rfweights("./worldclim.nc", nf, lon=lon_mat, lat=lat_mat, weightsfn="", fsmooth = TRUE)
 #' 
 #' # downscale using weights (ww=1. means do not use weights)
 #' r <- rainfarm(pr, lon_mat, lat_mat, nf, fsmooth = TRUE, fglob = FALSE, weights = ww, nens = 2, verbose = TRUE)
-- 
GitLab


From 9f27a8317b4f07358d74662cfb5ce95a5c40e05b Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Fri, 21 Dec 2018 09:06:36 +0100
Subject: [PATCH 07/48] further typos

---
 R/rainfarm.R | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/R/rainfarm.R b/R/rainfarm.R
index 3b9d2120..3aceaeb3 100644
--- a/R/rainfarm.R
+++ b/R/rainfarm.R
@@ -2,7 +2,7 @@
 #'
 #' @author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
 #'
-#' @description This function implements the RainFARM stochastic precipitation downscaling method.
+#' @description These functions implement the RainFARM stochastic precipitation downscaling method.
 #' Adapted for climate downscaling and including orographic correction.
 #' References: D'Onofrio et al. 2014, J of Hydrometeorology 15, 830-843; Rebora et. al 2006, JHM 7, 724
 #' The functions use julia_setup from the JuliaCall library, 
@@ -12,7 +12,6 @@
 #' @examples
 #' 
 #' source("rainfarm.R")
-#' source("rfweights.R")
 #' 
 #' nf <- 8   # Choose a downscaling by factor 8
 #' nens <- 3 # Number of ensemble members
-- 
GitLab


From c258b612a8b21ac535e31b6c8c1820b6212ac23b Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Fri, 21 Dec 2018 14:49:31 +0100
Subject: [PATCH 08/48] named list output

---
 R/rainfarm.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/rainfarm.R b/R/rainfarm.R
index 3aceaeb3..d3f1e68c 100644
--- a/R/rainfarm.R
+++ b/R/rainfarm.R
@@ -96,6 +96,6 @@ rainfarm <- function(pr, lon, lat, nf, fglob = FALSE,
                                 fglob = fglob, fsmooth = fsmooth,
                                 verbose = verbose, need_return = "R")))
   }
-
-  return(c(list(lon_f), list(lat_f), list(rd)))
+  return(list(lon = lon_f, lat = lat_f, data = rd))   
+# c(list(lon_f), list(lat_f), list(rd)))
 }
-- 
GitLab


From ffcbbb9d667ecbecfbeeb59e392f6040e8249a3c Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Fri, 21 Dec 2018 18:38:45 +0100
Subject: [PATCH 09/48] split rfweights and rainfarm, updated description

---
 R/rainfarm.R  | 61 ++++++++++++++-------------------------
 R/rfweights.R | 80 +++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 102 insertions(+), 39 deletions(-)
 create mode 100644 R/rfweights.R

diff --git a/R/rainfarm.R b/R/rainfarm.R
index d3f1e68c..28ab2fd9 100644
--- a/R/rainfarm.R
+++ b/R/rainfarm.R
@@ -2,17 +2,30 @@
 #'
 #' @author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
 #'
-#' @description These functions implement the RainFARM stochastic precipitation downscaling method.
+#' @description Perform RainFARM downscaling.
+#' These functions implement the RainFARM stochastic precipitation downscaling method.
 #' Adapted for climate downscaling and including orographic correction.
 #' References: D'Onofrio et al. 2014, J of Hydrometeorology 15, 830-843; Rebora et. al 2006, JHM 7, 724
 #' The functions use julia_setup from the JuliaCall library, 
 #' which takes a long time only the first time it is called on a machine
+#' 
+#' @param pr the precipitation array to downscal
+#' @param lon vector or array of longitudes
+#' @param lat vector or array of latitudes
+#' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
+#' @param slope prescribed or automatic spectral slope (default 0. == automatic)
+#' @param kmin first wavenumber for automatic spectral slope
+#' @param weights matrix (can be obtained using the \code{rfweights} function (default 1. == none)
+#' @param nens number of ensemble members to produce
+#' @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 verbose output of the Julia package
 #'
-#' @import JuliaCall, ncdf4
+#' @return A list containing the fine-scale longitudes, latitudes and the sequence of \code{nens} downscaled fields
+#'  
+#' @import JuliaCall
 #' @examples
 #' 
-#' source("rainfarm.R")
-#' 
 #' nf <- 8   # Choose a downscaling by factor 8
 #' nens <- 3 # Number of ensemble members
 #'
@@ -34,45 +47,16 @@
 #' 
 #' # downscale using weights (ww=1. means do not use weights)
 #' r <- rainfarm(pr, lon_mat, lat_mat, nf, fsmooth = TRUE, fglob = FALSE, weights = ww, nens = 2, verbose = TRUE)
-#' rd <- r[[1]]
-#' lon_f <- r[[2]]
-#' lat_f <- r[[3]]
-#'
+#' lon_f <- r$lon
+#' lat_f <- r$lat
+#' rd1 <- r$data[[1]]
+#' rd2 <- r$data[[2]]
 #' @export
 
-rfweights <- function(orofile, nf, lon=NULL, lat=NULL, reffile="",
-                      weightsfn="", varname="", fsmooth=TRUE) {
-
-  library(JuliaCall)
-  library(ncdf4)
-  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 <- "__reffile.nc"
-     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)
-
-  system("rm -f __reffile.nc gridrf_2_*.nc")
-  return(ww)
-}
-
 rainfarm <- function(pr, lon, lat, nf, fglob = FALSE,
                      fsmooth = TRUE, weights = 1., nens = 1,
                      verbose = FALSE, slope = 0, kmin=1) {
-  library(JuliaCall)
+#  library(JuliaCall)
   julia_setup(useRCall = FALSE)
   julia_library("RainFARM")
 
@@ -97,5 +81,4 @@ rainfarm <- function(pr, lon, lat, nf, fglob = FALSE,
                                 verbose = verbose, need_return = "R")))
   }
   return(list(lon = lon_f, lat = lat_f, data = rd))   
-# c(list(lon_f), list(lat_f), list(rd)))
 }
diff --git a/R/rfweights.R b/R/rfweights.R
new file mode 100644
index 00000000..bd3d341a
--- /dev/null
+++ b/R/rfweights.R
@@ -0,0 +1,80 @@
+#' RainFARM Stochastic Precipitation Downscaling
+#'
+#' @author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
+#'
+#' @description Compute orographic weights from a fine-scale precipitation climatology file.
+#' These functions implement the RainFARM stochastic precipitation downscaling method.
+#' Adapted for climate downscaling and including orographic correction.
+#' References: D'Onofrio et al. 2014, J of Hydrometeorology 15, 830-843; Rebora et. al 2006, JHM 7, 724
+#' The functions use julia_setup 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 
+#' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
+#' @param lon vector or array of longitudes (can be left out if \code{reffile} is used)
+#' @param lat vector or array of latitudes (can be left out 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 used
+#' @param varname optional name of variable to be read from reffile
+#' @param weightsfn optional filename to save the weights in netcdf format
+#' @param fsmooth logical to use smooth conservation (default) or large-scale box-average conservation
+#'
+#' @return The matrix containing the weights
+#'  
+#' @import JuliaCall, ncdf4
+#' @examples
+#'  
+#' 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(8,8,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
+#' 
+#' ww <- rfweights("./worldclim.nc", nf, lon=lon_mat, lat=lat_mat, weightsfn="", fsmooth = TRUE)
+#' 
+#' # downscale using weights (ww=1. means do not use weights)
+#' r <- rainfarm(pr, lon_mat, lat_mat, nf, fsmooth = TRUE, fglob = FALSE, weights = ww, nens = 2, verbose = TRUE)
+#' lon_f <- r$lon
+#' lat_f <- r$lat
+#' rd1 <- r$data[[1]]
+#' rd2 <- r$data[[2]]
+#' @export
+
+rfweights <- function(orofile, nf, lon=NULL, lat=NULL, reffile="",
+                      weightsfn="", varname="", fsmooth=TRUE) {
+
+#  library(JuliaCall)
+#  library(ncdf4)
+  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 <- "__reffile.nc"
+     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)
+
+  system("rm -f __reffile.nc gridrf_2_*.nc")
+  return(ww)
+}
-- 
GitLab


From 61264c246fe203d0445aeeacaaadc97273c1c31b Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Thu, 21 Feb 2019 11:50:43 +0200
Subject: [PATCH 10/48] first version with Nicolau's changes

---
 R/CST_RainFARM.R | 157 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 157 insertions(+)
 create mode 100644 R/CST_RainFARM.R

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
new file mode 100644
index 00000000..1fc36688
--- /dev/null
+++ b/R/CST_RainFARM.R
@@ -0,0 +1,157 @@
+#' RainFARM Stochastic Precipitation Downscaling
+#'
+#' @author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
+#'
+#' @description Perform RainFARM downscaling.
+#' These functions implement the RainFARM stochastic precipitation downscaling method.
+#' Adapted for climate downscaling and including orographic correction.
+#' References: D'Onofrio et al. 2014, J of Hydrometeorology 15, 830-843; Rebora et. al 2006, JHM 7, 724
+#' The functions use julia_setup from the JuliaCall library, 
+#' which takes a long time only the first time it is called on a machine
+#' 
+#' @param pr the precipitation array to downscal
+#' @param lon vector or array of longitudes
+#' @param lat vector or array of latitudes
+#' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
+#' @param slope prescribed or automatic spectral slope (default 0. == automatic)
+#' @param kmin first wavenumber for automatic spectral slope
+#' @param weights matrix (can be obtained using the \code{rfweights} function (default 1. == none)
+#' @param nens number of ensemble members to produce
+#' @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 verbose output of the Julia package
+#'
+#' @return A list containing the fine-scale longitudes, latitudes and the sequence of \code{nens} downscaled fields
+#'  
+#' @import JuliaCall
+#' @examples
+#' 
+#' 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(8,8,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
+#' 
+#' ww <- rfweights("./worldclim.nc", nf, lon=lon_mat, lat=lat_mat, weightsfn="", fsmooth = TRUE)
+#' 
+#' # downscale using weights (ww=1. means do not use weights)
+#' r <- rainfarm(pr, lon_mat, lat_mat, nf, fsmooth = TRUE, fglob = FALSE, weights = ww, nens = 2, verbose = TRUE)
+#' lon_f <- r$lon
+#' lat_f <- r$lat
+#' rd1 <- r$data[[1]]
+#' rd2 <- r$data[[2]]
+#' @export
+
+CST_RainFARM <- function(data, weights = 1., nf,
+                         slope = 0, kmin = 1,
+                         fglob = FALSE, fsmooth = FALSE,
+                         time_dim = "ftime",
+                         verbose = FALSE) {
+
+  res <- RainFARM(data$exp, data$lon, data$lat,
+           weights, nf, nens = 1,
+           slope, kmin,
+           fglob, fsmooth,
+           time_dim, drop_ens_dim = TRUE,
+           verbose)
+
+  n_extra_dims <- length(dim(res$data)) - 4
+  extra_dims <- NULL
+  if (n_extra_dims > 0) {
+    extra_dims <- 1:n_extra_dims + 3
+  }
+  data$exp <- multiApply:::.aperm2(res$data, c(3, 2, 1, extra_dims))
+
+  data$lon <- res$lon
+  data$lat <- res$lat
+  data$obs <- NULL
+
+  data
+}
+
+
+RainFARM <- function(data, lon, lat, weights = 1., nf,
+                     nens = 1, slope = 0, kmin = 1,
+                     fglob = FALSE, fsmooth = TRUE,
+                     time_dim = NULL, drop_ens_dim = FALSE,
+                     verbose = FALSE) {
+
+  # Check/detect time_dim
+  if (is.null(time_dim)) {
+    time_dim_names <- c("ftime", "sdate", "time")
+    time_dim <- which(time_dim_names %in% names(dim(data)))
+    if (length(time_dim) > 0) {
+      time_dim_name <- time_dim_names[time_dim[1]]
+    } else {
+      stop("Could not automatically detect a target time dimension ",
+           "in the provided data in 'data'.")
+    }
+  }
+
+  # 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]]
+
+  # Repeatedly apply .RainFARM
+  result <- Apply(data, c("lon", "lat", time_dim_name), .RainFARM,
+                  weights, nf, nens, slope, kmin,
+                  fglob, fsmooth, verbose)$output1
+
+  if (drop_ens_dim && nens == 1) {
+    dim(result) <- dim(result)[-which(names(dim(result)) == "member")[1]]
+  } else {
+    if (drop_ens_dim) {
+      stop("Requested 'drop_ens_dim = TRUE' but 'nens > 1'. ",
+           "Preserving 'member' dimension.")
+    }
+  }
+
+  list(data = result, lon = lon_f, lat = lat_f)
+}
+
+
+.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(lon = dim(weights)[1],
+                   lat = dim(weights)[2],
+                   dim(pr)[3],
+                   member = nens)
+  names(result_dims)[3] <- names(dim(pr))[3]
+  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")
+  }
+}
-- 
GitLab


From 6e2ee26b761728fd5a67f4daedffcf554aff2568 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Mon, 4 Mar 2019 00:34:58 +0100
Subject: [PATCH 11/48] rewritten for CSTools compliance according to Nicolau's
 template

---
 R/CST_RainFARM.R | 23 ++++++-------
 R/rainfarm.R     | 84 ------------------------------------------------
 2 files changed, 12 insertions(+), 95 deletions(-)
 delete mode 100644 R/rainfarm.R

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index 1fc36688..b52ab7ce 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -54,19 +54,18 @@
 #' @export
 
 CST_RainFARM <- function(data, weights = 1., nf,
-                         slope = 0, kmin = 1,
+                         slope = 0, kmin = 1, nens = 1,
                          fglob = FALSE, fsmooth = FALSE,
-                         time_dim = "ftime",
+                         time_dim = NULL,
                          verbose = FALSE) {
 
   res <- RainFARM(data$exp, data$lon, data$lat,
-           weights, nf, nens = 1,
+           weights, nf, nens,
            slope, kmin,
            fglob, fsmooth,
            time_dim, drop_ens_dim = TRUE,
            verbose)
-
-  n_extra_dims <- length(dim(res$data)) - 4
+  n_extra_dims <- length(dim(res$data)) - 3
   extra_dims <- NULL
   if (n_extra_dims > 0) {
     extra_dims <- 1:n_extra_dims + 3
@@ -77,7 +76,7 @@ CST_RainFARM <- function(data, weights = 1., nf,
   data$lat <- res$lat
   data$obs <- NULL
 
-  data
+  return(data)
 }
 
 
@@ -90,9 +89,9 @@ RainFARM <- function(data, lon, lat, weights = 1., nf,
   # Check/detect time_dim
   if (is.null(time_dim)) {
     time_dim_names <- c("ftime", "sdate", "time")
-    time_dim <- which(time_dim_names %in% names(dim(data)))
-    if (length(time_dim) > 0) {
-      time_dim_name <- time_dim_names[time_dim[1]]
+    time_dim_num <- which(time_dim_names %in% names(dim(data)))
+    if (length(time_dim_num) > 0) {
+      time_dim <- time_dim_names[time_dim_num[1]]
     } else {
       stop("Could not automatically detect a target time dimension ",
            "in the provided data in 'data'.")
@@ -115,7 +114,7 @@ RainFARM <- function(data, lon, lat, weights = 1., nf,
   lat_f <- r[[2]]
 
   # Repeatedly apply .RainFARM
-  result <- Apply(data, c("lon", "lat", time_dim_name), .RainFARM,
+  result <- Apply(data, c("lon", "lat", time_dim), .RainFARM,
                   weights, nf, nens, slope, kmin,
                   fglob, fsmooth, verbose)$output1
 
@@ -128,12 +127,13 @@ RainFARM <- function(data, lon, lat, weights = 1., nf,
     }
   }
 
-  list(data = result, lon = lon_f, lat = lat_f)
+  return(list(data = result, lon = lon_f, lat = lat_f))
 }
 
 
 .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),
@@ -154,4 +154,5 @@ RainFARM <- function(data, lon, lat, weights = 1., nf,
                              fglob = fglob, fsmooth = fsmooth,
                              verbose = verbose, need_return = "R")
   }
+  return(r)
 }
diff --git a/R/rainfarm.R b/R/rainfarm.R
deleted file mode 100644
index 28ab2fd9..00000000
--- a/R/rainfarm.R
+++ /dev/null
@@ -1,84 +0,0 @@
-#' RainFARM Stochastic Precipitation Downscaling
-#'
-#' @author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
-#'
-#' @description Perform RainFARM downscaling.
-#' These functions implement the RainFARM stochastic precipitation downscaling method.
-#' Adapted for climate downscaling and including orographic correction.
-#' References: D'Onofrio et al. 2014, J of Hydrometeorology 15, 830-843; Rebora et. al 2006, JHM 7, 724
-#' The functions use julia_setup from the JuliaCall library, 
-#' which takes a long time only the first time it is called on a machine
-#' 
-#' @param pr the precipitation array to downscal
-#' @param lon vector or array of longitudes
-#' @param lat vector or array of latitudes
-#' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
-#' @param slope prescribed or automatic spectral slope (default 0. == automatic)
-#' @param kmin first wavenumber for automatic spectral slope
-#' @param weights matrix (can be obtained using the \code{rfweights} function (default 1. == none)
-#' @param nens number of ensemble members to produce
-#' @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 verbose output of the Julia package
-#'
-#' @return A list containing the fine-scale longitudes, latitudes and the sequence of \code{nens} downscaled fields
-#'  
-#' @import JuliaCall
-#' @examples
-#' 
-#' 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(8,8,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
-#' 
-#' ww <- rfweights("./worldclim.nc", nf, lon=lon_mat, lat=lat_mat, weightsfn="", fsmooth = TRUE)
-#' 
-#' # downscale using weights (ww=1. means do not use weights)
-#' r <- rainfarm(pr, lon_mat, lat_mat, nf, fsmooth = TRUE, fglob = FALSE, weights = ww, nens = 2, verbose = TRUE)
-#' lon_f <- r$lon
-#' lat_f <- r$lat
-#' rd1 <- r$data[[1]]
-#' rd2 <- r$data[[2]]
-#' @export
-
-rainfarm <- function(pr, lon, lat, nf, fglob = FALSE,
-                     fsmooth = TRUE, weights = 1., nens = 1,
-                     verbose = FALSE, slope = 0, kmin=1) {
-#  library(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]]
-
-  if (slope == 0) {
-    r <- julia_call("fft3d", pr, need_return = "R")
-    fxp <- r[[1]]
-    sx <- julia_call("fitslopex", fxp, kmin = as.integer(kmin),
-                     need_return = "R")
-  } else {
-    sx <- slope
-  }
-
-  rd <- list()
-  for (j in 1:nens) {
-    rd <- c(rd, list(julia_call("rainfarm", pr, sx, as.integer(nf), weights,
-                                fglob = fglob, fsmooth = fsmooth,
-                                verbose = verbose, need_return = "R")))
-  }
-  return(list(lon = lon_f, lat = lat_f, data = rd))   
-}
-- 
GitLab


From d6d971cca3088959ed2898707c34f81277f94455 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Thu, 14 Mar 2019 22:49:31 +0100
Subject: [PATCH 12/48] upgraded documentation

---
 R/CST_RainFARM.R | 118 +++++++++++++++++++++++++++++++----------------
 R/rfweights.R    |  49 ++++++++------------
 2 files changed, 95 insertions(+), 72 deletions(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index b52ab7ce..61bd7b66 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -1,4 +1,6 @@
-#' RainFARM Stochastic Precipitation Downscaling
+#' @name rainfarm
+#' @rdname rainfarm
+#' @title RainFARM Stochastic Precipitation Downscaling
 #'
 #' @author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
 #'
@@ -8,77 +10,107 @@
 #' References: D'Onofrio et al. 2014, J of Hydrometeorology 15, 830-843; Rebora et. al 2006, JHM 7, 724
 #' The functions use julia_setup from the JuliaCall library, 
 #' which takes a long time only the first time it is called on a machine
-#' 
-#' @param pr the precipitation array to downscal
-#' @param lon vector or array of longitudes
-#' @param lat vector or array of latitudes
+#'
+#' @param weights matrix (can be obtained using the \code{rfweights} function (default 1. == none)
 #' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
 #' @param slope prescribed or automatic spectral slope (default 0. == automatic)
 #' @param kmin first wavenumber for automatic spectral slope
-#' @param weights matrix (can be obtained using the \code{rfweights} function (default 1. == none)
 #' @param nens number of ensemble members to produce
 #' @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 verbose output of the Julia package
-#'
-#' @return A list containing the fine-scale longitudes, latitudes and the sequence of \code{nens} downscaled fields
+#' @param time_dim name of time dimension (e.g. "ftime", "sdate", "time" ...)
+#' @param verbose logical for verbose output of the Julia package
 #'  
 #' @import JuliaCall
+NULL
+
+#' RainFARM Stochastic Precipitation Downscaling for CSTools objects
+#' @rdname rainfarm
+#' @param s2d_data S2D list containing the precipitation array to downscale
+#' @return CST_RainFARM() returns a downscaled S2D object
+#' @export
 #' @examples
-#' 
+#' #Example using CST_RainFARM for a CSTools object
 #' 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(8,8,20)
-#' lon_mat <- seq(10,13.5,0.5) # could also be a 2d matrix
-#' lat_mat <- seq(40,43.5,0.5)
-#'
+#' 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))
+#' s2dv_data <- list(exp = exp, lon = lon, lat = lat)
 #' # 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
-#' 
-#' ww <- rfweights("./worldclim.nc", nf, lon=lon_mat, lat=lat_mat, weightsfn="", fsmooth = TRUE)
-#' 
-#' # downscale using weights (ww=1. means do not use weights)
-#' r <- rainfarm(pr, lon_mat, lat_mat, nf, fsmooth = TRUE, fglob = FALSE, weights = ww, nens = 2, verbose = TRUE)
-#' lon_f <- r$lon
-#' lat_f <- r$lat
-#' rd1 <- r$data[[1]]
-#' rd2 <- r$data[[2]]
-#' @export
+#' res_s2dv <- CST_RainFARM(s2dv_data, ww, nf)
+#' str(res_s2dv)
+#' #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_s2dv$exp)
+#' #  ftime     lat     lon dataset  member   sdate 
+#' #      4      64      64       1       2       3
+#'
 
-CST_RainFARM <- function(data, weights = 1., nf,
+CST_RainFARM <- function(s2d_data, weights = 1., nf,
                          slope = 0, kmin = 1, nens = 1,
                          fglob = FALSE, fsmooth = FALSE,
                          time_dim = NULL,
                          verbose = FALSE) {
 
-  res <- RainFARM(data$exp, data$lon, data$lat,
+  res <- RainFARM(s2d_data$exp, s2d_data$lon, s2d_data$lat,
            weights, nf, nens,
            slope, kmin,
            fglob, fsmooth,
            time_dim, drop_ens_dim = TRUE,
            verbose)
-  n_extra_dims <- length(dim(res$data)) - 3
+  n_extra_dims <- length(dim(res$s2d_data)) - 3
   extra_dims <- NULL
   if (n_extra_dims > 0) {
     extra_dims <- 1:n_extra_dims + 3
   }
-  data$exp <- multiApply:::.aperm2(res$data, c(3, 2, 1, extra_dims))
+  s2d_data$exp <- multiApply:::.aperm2(res$s2d_data, c(3, 2, 1, extra_dims))
 
-  data$lon <- res$lon
-  data$lat <- res$lat
-  data$obs <- NULL
+  s2d_data$lon <- res$lon
+  s2d_data$lat <- res$lat
+  s2d_data$obs <- NULL
 
-  return(data)
+  return(s2d_data)
 }
 
+#' Reduced RainFARM function
+#' @rdname rainfarm
+#' @param data precipitation array to downscale
+#' @param lon vector or array of longitudes
+#' @param lat vector or array of latitudes
+#' @param drop_ens_dim logical to remove the ensemble diemnsion if nens==1
+#' @return RainFARM() returns a list containing the fine-scale longitudes, latitudes and the sequence of \code{nens} downscaled fields
+#' @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
+#' ww <- rfweights("./worldclim.nc", nf, lon=lon_mat, lat=lat_mat, weightsfn="", fsmooth = TRUE)
+#' # downscale using weights (ww=1. means do not use weights)
+#' r <- 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)
+#' #member  ftime	lat	lon
+#' # 	3 	20 	64 	64
 
 RainFARM <- function(data, lon, lat, weights = 1., nf,
                      nens = 1, slope = 0, kmin = 1,
@@ -130,6 +162,10 @@ RainFARM <- function(data, lon, lat, weights = 1., nf,
   return(list(data = result, lon = lon_f, lat = lat_f))
 }
 
+#' Atomic RainFARM
+#' @param pr precipitation array to downscale
+#' @return .RainFARM returns a downscaled array
+#' @noRd
 
 .RainFARM <- function(pr, weights, nf, nens, slope, kmin,
                       fglob, fsmooth, verbose) {
diff --git a/R/rfweights.R b/R/rfweights.R
index bd3d341a..7c1ab28b 100644
--- a/R/rfweights.R
+++ b/R/rfweights.R
@@ -1,12 +1,12 @@
-#' RainFARM Stochastic Precipitation Downscaling
+#' RainFARM Stochastic Precipitation Downscaling Weights
 #'
 #' @author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
 #'
 #' @description Compute orographic weights from a fine-scale precipitation climatology file.
-#' These functions implement the RainFARM stochastic precipitation downscaling method.
-#' Adapted for climate downscaling and including orographic correction.
-#' References: D'Onofrio et al. 2014, J of Hydrometeorology 15, 830-843; Rebora et. al 2006, JHM 7, 724
-#' The functions use julia_setup from the JuliaCall library, 
+#'
+#' 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
+#'
+#' The rfweights, CST_RainFARM, RainFARM functions use julia_setup 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 
@@ -18,36 +18,23 @@
 #' @param weightsfn optional filename to save the weights in netcdf format
 #' @param fsmooth logical to use smooth conservation (default) or large-scale box-average conservation
 #'
-#' @return The matrix containing the weights
+#' @return The matrix containing the weights with dimensions (lon, lat)
 #'  
 #' @import JuliaCall, ncdf4
 #' @examples
 #'  
-#' 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(8,8,20)
+#' # Create weights to be used with CST_RainFARM() or RainFARM()
+#' # Specify a weightsfn filename if you wish to save the weights
+#' # Using an external fine-scale climatology file.
+#' 
+#' # 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)
-#'
-#' # 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
-#' 
+#' nf <- 8
 #' ww <- rfweights("./worldclim.nc", nf, lon=lon_mat, lat=lat_mat, weightsfn="", fsmooth = TRUE)
-#' 
-#' # downscale using weights (ww=1. means do not use weights)
-#' r <- rainfarm(pr, lon_mat, lat_mat, nf, fsmooth = TRUE, fglob = FALSE, weights = ww, nens = 2, verbose = TRUE)
-#' lon_f <- r$lon
-#' lat_f <- r$lat
-#' rd1 <- r$data[[1]]
-#' rd2 <- r$data[[2]]
+#'
+#' # or use a reference file (e.g. the input file itself) to specify the coordinates
+#' ww <- rfweights("./worldclim.nc", nf, reffile="./input.nc", weightsfn="", fsmooth = TRUE)
 #' @export
 
 rfweights <- function(orofile, nf, lon=NULL, lat=NULL, reffile="",
@@ -67,14 +54,14 @@ rfweights <- function(orofile, nf, lon=NULL, lat=NULL, reffile="",
                        unlim = T)
      my_ncdf <- ncvar_def("grid", "m", list(xvar, yvar, tvar), -999,
                           longname = "grid", prec = "single")
-     reffile <- "__reffile.nc"
+     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)
-
-  system("rm -f __reffile.nc gridrf_2_*.nc")
+  
+  unlink(reffile)
   return(ww)
 }
-- 
GitLab


From 9c39b9383f6297727961af3a84b919e8cd85ce8a Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Thu, 14 Mar 2019 23:02:52 +0100
Subject: [PATCH 13/48] remove extra temporary file which julia RF does not
 remove

---
 R/rfweights.R | 1 +
 1 file changed, 1 insertion(+)

diff --git a/R/rfweights.R b/R/rfweights.R
index 7c1ab28b..27a31bcf 100644
--- a/R/rfweights.R
+++ b/R/rfweights.R
@@ -63,5 +63,6 @@ rfweights <- function(orofile, nf, lon=NULL, lat=NULL, reffile="",
                    weightsfn = weightsfn, varname = varname, fsmooth = fsmooth)
   
   unlink(reffile)
+  system("rm -f gridrf_2_*.nc")
   return(ww)
 }
-- 
GitLab


From 171c7e6d26ffb66f0807ef48612244d8bf731ed7 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Thu, 14 Mar 2019 23:19:47 +0100
Subject: [PATCH 14/48] typo in imports

---
 R/rfweights.R | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/R/rfweights.R b/R/rfweights.R
index 27a31bcf..81699a7b 100644
--- a/R/rfweights.R
+++ b/R/rfweights.R
@@ -20,7 +20,8 @@
 #'
 #' @return The matrix containing the weights with dimensions (lon, lat)
 #'  
-#' @import JuliaCall, ncdf4
+#' @import ncdf4
+#' @import JuliaCall
 #' @examples
 #'  
 #' # Create weights to be used with CST_RainFARM() or RainFARM()
-- 
GitLab


From 6b6b57415728438db07223a23caac2039705b21f Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Thu, 14 Mar 2019 23:44:52 +0100
Subject: [PATCH 15/48] removal of temp file not needed since Julia RainFARM
 has been fixed

---
 R/rfweights.R | 1 -
 1 file changed, 1 deletion(-)

diff --git a/R/rfweights.R b/R/rfweights.R
index 81699a7b..68aaf104 100644
--- a/R/rfweights.R
+++ b/R/rfweights.R
@@ -64,6 +64,5 @@ rfweights <- function(orofile, nf, lon=NULL, lat=NULL, reffile="",
                    weightsfn = weightsfn, varname = varname, fsmooth = fsmooth)
   
   unlink(reffile)
-  system("rm -f gridrf_2_*.nc")
   return(ww)
 }
-- 
GitLab


From 3f15cb90e9bd3d21620b3ab6076d70e67bb6724e Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Fri, 15 Mar 2019 15:21:55 +0100
Subject: [PATCH 16/48] reverted bug introduced in doc update

---
 R/CST_RainFARM.R | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index 61bd7b66..60073681 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -57,19 +57,18 @@ CST_RainFARM <- function(s2d_data, weights = 1., nf,
                          fglob = FALSE, fsmooth = FALSE,
                          time_dim = NULL,
                          verbose = FALSE) {
-
   res <- RainFARM(s2d_data$exp, s2d_data$lon, s2d_data$lat,
            weights, nf, nens,
            slope, kmin,
            fglob, fsmooth,
            time_dim, drop_ens_dim = TRUE,
            verbose)
-  n_extra_dims <- length(dim(res$s2d_data)) - 3
+  n_extra_dims <- length(dim(res$data)) - 3
   extra_dims <- NULL
   if (n_extra_dims > 0) {
     extra_dims <- 1:n_extra_dims + 3
   }
-  s2d_data$exp <- multiApply:::.aperm2(res$s2d_data, c(3, 2, 1, extra_dims))
+  s2d_data$exp <- multiApply:::.aperm2(res$data, c(3, 2, 1, extra_dims))
 
   s2d_data$lon <- res$lon
   s2d_data$lat <- res$lat
-- 
GitLab


From d7c032d7e2ebba0a925ba79ea98492b73383444e Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Fri, 15 Mar 2019 15:50:01 +0100
Subject: [PATCH 17/48] further improvement to docs

---
 R/CST_RainFARM.R | 52 ++++++++++++++++++++++++++++++++----------------
 R/rfweights.R    |  8 +++-----
 2 files changed, 38 insertions(+), 22 deletions(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index 60073681..b9d328cf 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -1,16 +1,15 @@
-#' @name rainfarm
-#' @rdname rainfarm
-#' @title RainFARM Stochastic Precipitation Downscaling
+#' @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 Perform RainFARM downscaling.
-#' These functions implement the RainFARM stochastic precipitation downscaling method.
+#' @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.
-#' References: D'Onofrio et al. 2014, J of Hydrometeorology 15, 830-843; Rebora et. al 2006, JHM 7, 724
-#' The functions use julia_setup from the JuliaCall library, 
-#' which takes a long time only the first time it is called on a machine
+#' 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 julia_setup from the JuliaCall library, 
+#' which takes a long time only the first time it is called on a machine.
 #'
+#' @param s2d_data CSTools object containing the precipitation array to downscale
 #' @param weights matrix (can be obtained using the \code{rfweights} function (default 1. == none)
 #' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
 #' @param slope prescribed or automatic spectral slope (default 0. == automatic)
@@ -22,12 +21,7 @@
 #' @param verbose logical for verbose output of the Julia package
 #'  
 #' @import JuliaCall
-NULL
-
-#' RainFARM Stochastic Precipitation Downscaling for CSTools objects
-#' @rdname rainfarm
-#' @param s2d_data S2D list containing the precipitation array to downscale
-#' @return CST_RainFARM() returns a downscaled S2D object
+#' @return CST_RainFARM() returns a downscaled CSTools object.
 #' @export
 #' @examples
 #' #Example using CST_RainFARM for a CSTools object
@@ -77,13 +71,29 @@ CST_RainFARM <- function(s2d_data, weights = 1., nf,
   return(s2d_data)
 }
 
-#' Reduced RainFARM function
-#' @rdname rainfarm
+#' @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 an array with dims (lon, lat, time) in input.
+#' 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 julia_setup 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
 #' @param lon vector or array of longitudes
 #' @param lat vector or array of latitudes
+#' @param weights matrix (can be obtained using the \code{rfweights} function (default 1. == none)
+#' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
+#' @param slope prescribed or automatic spectral slope (default 0. == automatic)
+#' @param kmin first wavenumber for automatic spectral slope
+#' @param nens number of ensemble members to produce
+#' @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 time_dim name of time dimension (e.g. "ftime", "sdate", "time" ...)
+#' @param verbose logical for verbose output of the Julia package
 #' @param drop_ens_dim logical to remove the ensemble diemnsion if nens==1
-#' @return RainFARM() returns a list containing the fine-scale longitudes, latitudes and the sequence of \code{nens} downscaled fields
+#' @return RainFARM() returns a list containing the fine-scale longitudes, latitudes and the sequence of \code{nens} downscaled fields.
+#' @import JuliaCall
 #' @export
 #' @examples
 #' # Example for the 'reduced' RainFARM function 
@@ -163,6 +173,14 @@ RainFARM <- function(data, lon, lat, weights = 1., nf,
 
 #' Atomic RainFARM
 #' @param pr precipitation array to downscale
+#' @param weights matrix (can be obtained using the \code{rfweights} function (default 1. == none)
+#' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
+#' @param slope prescribed or automatic spectral slope (default 0. == automatic)
+#' @param kmin first wavenumber for automatic spectral slope
+#' @param nens number of ensemble members to produce
+#' @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
 #' @return .RainFARM returns a downscaled array
 #' @noRd
 
diff --git a/R/rfweights.R b/R/rfweights.R
index 68aaf104..5cf841ea 100644
--- a/R/rfweights.R
+++ b/R/rfweights.R
@@ -4,9 +4,9 @@
 #'
 #' @description Compute 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
+#' 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 .
 #'
-#' The rfweights, CST_RainFARM, RainFARM functions use julia_setup from the JuliaCall library, 
+#' Please notice: The rfweights() function uses julia_setup 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 
@@ -17,9 +17,7 @@
 #' @param varname optional name of variable to be read from reffile
 #' @param weightsfn optional filename to save the weights in netcdf format
 #' @param fsmooth logical to use smooth conservation (default) or large-scale box-average conservation
-#'
-#' @return The matrix containing the weights with dimensions (lon, lat)
-#'  
+#' @return The matrix containing the weights with dimensions (lon, lat).
 #' @import ncdf4
 #' @import JuliaCall
 #' @examples
-- 
GitLab


From 0fdbd752a69ca81360391d407cb6a30e7668e45c Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Fri, 15 Mar 2019 16:20:43 +0100
Subject: [PATCH 18/48] make lintr happy

---
 R/rfweights.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/rfweights.R b/R/rfweights.R
index 5cf841ea..0c104dd8 100644
--- a/R/rfweights.R
+++ b/R/rfweights.R
@@ -60,7 +60,7 @@ rfweights <- function(orofile, nf, lon=NULL, lat=NULL, reffile="",
 
   ww <- julia_call("rfweights", orofile, reffile, as.integer(nf),
                    weightsfn = weightsfn, varname = varname, fsmooth = fsmooth)
-  
+
   unlink(reffile)
   return(ww)
 }
-- 
GitLab


From db934e092a2aac41557dea3d2c82a278b6ffe0c4 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Fri, 15 Mar 2019 16:46:25 +0100
Subject: [PATCH 19/48] renamed rfweights to RFWeights

---
 R/CST_RainFARM.R               | 18 +++++++++---------
 R/{rfweights.R => RFWeights.R} |  8 ++++----
 2 files changed, 13 insertions(+), 13 deletions(-)
 rename R/{rfweights.R => RFWeights.R} (92%)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index b9d328cf..fbab5d70 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -10,7 +10,7 @@
 #' which takes a long time only the first time it is called on a machine.
 #'
 #' @param s2d_data CSTools object containing the precipitation array to downscale
-#' @param weights matrix (can be obtained using the \code{rfweights} function (default 1. == none)
+#' @param weights matrix (can be obtained using the \code{RFWeights} function (default 1. == none)
 #' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
 #' @param slope prescribed or automatic spectral slope (default 0. == automatic)
 #' @param kmin first wavenumber for automatic spectral slope
@@ -35,7 +35,7 @@
 #' s2dv_data <- list(exp = exp, lon = lon, lat = lat)
 #' # Create a test array of weights
 #' ww <- array(1., dim = c(8 * nf, 8 * nf))
-#' res_s2dv <- CST_RainFARM(s2dv_data, ww, nf)
+#' res_s2dv <- CST_RainFARM(s2dv_data, nf, ww)
 #' str(res_s2dv)
 #' #List of 3
 #' # $ exp: num [1:4, 1:64, 1:64, 1, 1:2, 1:3] 201 115 119 307 146 ...
@@ -46,13 +46,13 @@
 #' #      4      64      64       1       2       3
 #'
 
-CST_RainFARM <- function(s2d_data, weights = 1., nf,
+CST_RainFARM <- function(s2d_data, nf, weights = 1.,
                          slope = 0, kmin = 1, nens = 1,
                          fglob = FALSE, fsmooth = FALSE,
                          time_dim = NULL,
                          verbose = FALSE) {
   res <- RainFARM(s2d_data$exp, s2d_data$lon, s2d_data$lat,
-           weights, nf, nens,
+           nf, weights, nens,
            slope, kmin,
            fglob, fsmooth,
            time_dim, drop_ens_dim = TRUE,
@@ -82,7 +82,7 @@ CST_RainFARM <- function(s2d_data, weights = 1., nf,
 #' @param data precipitation array to downscale
 #' @param lon vector or array of longitudes
 #' @param lat vector or array of latitudes
-#' @param weights matrix (can be obtained using the \code{rfweights} function (default 1. == none)
+#' @param weights matrix (can be obtained using the \code{RFWeights} function (default 1. == none)
 #' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
 #' @param slope prescribed or automatic spectral slope (default 0. == automatic)
 #' @param kmin first wavenumber for automatic spectral slope
@@ -109,9 +109,9 @@ CST_RainFARM <- function(s2d_data, weights = 1., nf,
 #' 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)
+#' ww <- RFWeights("./worldclim.nc", nf, lon=lon_mat, lat=lat_mat, weightsfn="", fsmooth = TRUE)
 #' # downscale using weights (ww=1. means do not use weights)
-#' r <- rainfarm(pr, lon_mat, lat_mat, nf, fsmooth = TRUE, fglob = FALSE, weights = ww, nens = 2, verbose = TRUE)
+#' r <- 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 ...
@@ -121,7 +121,7 @@ CST_RainFARM <- function(s2d_data, weights = 1., nf,
 #' #member  ftime	lat	lon
 #' # 	3 	20 	64 	64
 
-RainFARM <- function(data, lon, lat, weights = 1., nf,
+RainFARM <- function(data, lon, lat, nf, weights = 1.,
                      nens = 1, slope = 0, kmin = 1,
                      fglob = FALSE, fsmooth = TRUE,
                      time_dim = NULL, drop_ens_dim = FALSE,
@@ -173,7 +173,7 @@ RainFARM <- function(data, lon, lat, weights = 1., nf,
 
 #' Atomic RainFARM
 #' @param pr precipitation array to downscale
-#' @param weights matrix (can be obtained using the \code{rfweights} function (default 1. == none)
+#' @param weights matrix (can be obtained using the \code{RFWeights} function (default 1. == none)
 #' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
 #' @param slope prescribed or automatic spectral slope (default 0. == automatic)
 #' @param kmin first wavenumber for automatic spectral slope
diff --git a/R/rfweights.R b/R/RFWeights.R
similarity index 92%
rename from R/rfweights.R
rename to R/RFWeights.R
index 0c104dd8..d6eab6e1 100644
--- a/R/rfweights.R
+++ b/R/RFWeights.R
@@ -6,7 +6,7 @@
 #'
 #' 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 rfweights() function uses julia_setup from the JuliaCall library, 
+#' Please notice: The RFWeights() function uses julia_setup 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 
@@ -30,13 +30,13 @@
 #' 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)
+#' 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", weightsfn="", fsmooth = TRUE)
+#' ww <- RFWeights("./worldclim.nc", nf, reffile="./input.nc", weightsfn="", fsmooth = TRUE)
 #' @export
 
-rfweights <- function(orofile, nf, lon=NULL, lat=NULL, reffile="",
+RFWeights <- function(orofile, nf, lon=NULL, lat=NULL, reffile="",
                       weightsfn="", varname="", fsmooth=TRUE) {
 
 #  library(JuliaCall)
-- 
GitLab


From 047dd9b016523a77ec8f9fc04c5b101a2879eab6 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Fri, 15 Mar 2019 17:04:33 +0100
Subject: [PATCH 20/48] added import of multiApply

---
 R/CST_RainFARM.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index fbab5d70..8d50f17e 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -20,7 +20,6 @@
 #' @param time_dim name of time dimension (e.g. "ftime", "sdate", "time" ...)
 #' @param verbose logical for verbose output of the Julia package
 #'  
-#' @import JuliaCall
 #' @return CST_RainFARM() returns a downscaled CSTools object.
 #' @export
 #' @examples
@@ -94,6 +93,7 @@ CST_RainFARM <- function(s2d_data, nf, weights = 1.,
 #' @param drop_ens_dim logical to remove the ensemble diemnsion if nens==1
 #' @return RainFARM() returns a list containing the fine-scale longitudes, latitudes and the sequence of \code{nens} downscaled fields.
 #' @import JuliaCall
+#' @import multiApply
 #' @export
 #' @examples
 #' # Example for the 'reduced' RainFARM function 
-- 
GitLab


From ef713aced9e0ebcfab7c99ac3754389f430243b9 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Fri, 15 Mar 2019 17:22:29 +0100
Subject: [PATCH 21/48] bug when weights scalar

---
 R/CST_RainFARM.R | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index 8d50f17e..02473444 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -21,6 +21,7 @@
 #' @param verbose logical for verbose output of the Julia package
 #'  
 #' @return CST_RainFARM() returns a downscaled CSTools object.
+#' @import multiApply
 #' @export
 #' @examples
 #' #Example using CST_RainFARM for a CSTools object
@@ -195,8 +196,8 @@ RainFARM <- function(data, lon, lat, nf, weights = 1.,
     sx <- slope
   }
 
-  result_dims <- c(lon = dim(weights)[1],
-                   lat = dim(weights)[2],
+  result_dims <- c(lon = dim(pr)[1] * nf,
+                   lat = dim(pr)[2] * nf,
                    dim(pr)[3],
                    member = nens)
   names(result_dims)[3] <- names(dim(pr))[3]
-- 
GitLab


From 610210698bc879659bf5f2e4ee6c975bbb191756 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Fri, 15 Mar 2019 23:32:35 +0100
Subject: [PATCH 22/48] added named dimensions to RFWeights

---
 R/RFWeights.R | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/R/RFWeights.R b/R/RFWeights.R
index d6eab6e1..1d2de5a6 100644
--- a/R/RFWeights.R
+++ b/R/RFWeights.R
@@ -39,8 +39,6 @@
 RFWeights <- function(orofile, nf, lon=NULL, lat=NULL, reffile="",
                       weightsfn="", varname="", fsmooth=TRUE) {
 
-#  library(JuliaCall)
-#  library(ncdf4)
   julia_setup(useRCall = FALSE)
   julia_library("RainFARM")
   if ( (reffile == "") & ( (typeof(lon) == "NULL") | (typeof(lat) == "NULL"))) {
@@ -61,6 +59,7 @@ RFWeights <- function(orofile, nf, lon=NULL, lat=NULL, reffile="",
   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)
 }
-- 
GitLab


From d37b614570f4d229e5e488447a908c347069237f Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Mon, 18 Mar 2019 22:05:20 +0100
Subject: [PATCH 23/48] drop_ens_dim fix

---
 R/CST_RainFARM.R | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index 02473444..e44302a2 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -19,6 +19,7 @@
 #' @param fsmooth logical to conserve precipitation with a smnoothing kernel (default TRUE)
 #' @param time_dim name of time dimension (e.g. "ftime", "sdate", "time" ...)
 #' @param verbose logical for verbose output of the Julia package
+#' @param drop_ens_dim logical to remove the ensemble dimension if nens==1 (default FALSE)
 #'  
 #' @return CST_RainFARM() returns a downscaled CSTools object.
 #' @import multiApply
@@ -49,13 +50,14 @@
 CST_RainFARM <- function(s2d_data, nf, weights = 1.,
                          slope = 0, kmin = 1, nens = 1,
                          fglob = FALSE, fsmooth = FALSE,
-                         time_dim = NULL,
-                         verbose = FALSE) {
+                         time_dim = NULL, verbose = FALSE,
+                         drop_ens_dim = FALSE) {
+
   res <- RainFARM(s2d_data$exp, s2d_data$lon, s2d_data$lat,
            nf, weights, nens,
            slope, kmin,
            fglob, fsmooth,
-           time_dim, drop_ens_dim = TRUE,
+           time_dim, drop_ens_dim,
            verbose)
   n_extra_dims <- length(dim(res$data)) - 3
   extra_dims <- NULL
@@ -91,7 +93,7 @@ CST_RainFARM <- function(s2d_data, nf, weights = 1.,
 #' @param fsmooth logical to conserve precipitation with a smnoothing kernel (default TRUE)
 #' @param time_dim name of time dimension (e.g. "ftime", "sdate", "time" ...)
 #' @param verbose logical for verbose output of the Julia package
-#' @param drop_ens_dim logical to remove the ensemble diemnsion if nens==1
+#' @param drop_ens_dim logical to remove the ensemble dimension if nens==1
 #' @return RainFARM() returns a list containing the fine-scale longitudes, latitudes and the sequence of \code{nens} downscaled fields.
 #' @import JuliaCall
 #' @import multiApply
-- 
GitLab


From b5f03f9bb56afef613152974aff714d197ad2d3c Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Mon, 18 Mar 2019 22:25:36 +0100
Subject: [PATCH 24/48] ignore drop_ens_dim if nens>1

---
 R/CST_RainFARM.R | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index e44302a2..95218bd3 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -164,11 +164,6 @@ RainFARM <- function(data, lon, lat, nf, weights = 1.,
 
   if (drop_ens_dim && nens == 1) {
     dim(result) <- dim(result)[-which(names(dim(result)) == "member")[1]]
-  } else {
-    if (drop_ens_dim) {
-      stop("Requested 'drop_ens_dim = TRUE' but 'nens > 1'. ",
-           "Preserving 'member' dimension.")
-    }
   }
 
   return(list(data = result, lon = lon_f, lat = lat_f))
-- 
GitLab


From 9ee53ab3ddadb503aec725990e8c11da651f1beb Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Tue, 19 Mar 2019 00:18:57 +0100
Subject: [PATCH 25/48] fix lon and lat names

---
 R/CST_RainFARM.R | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index 95218bd3..d45c6ff4 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -158,6 +158,7 @@ RainFARM <- function(data, lon, lat, nf, weights = 1.,
   lat_f <- r[[2]]
 
   # Repeatedly apply .RainFARM
+  print(dim(data))
   result <- Apply(data, c("lon", "lat", time_dim), .RainFARM,
                   weights, nf, nens, slope, kmin,
                   fglob, fsmooth, verbose)$output1
@@ -193,11 +194,7 @@ RainFARM <- function(data, lon, lat, nf, weights = 1.,
     sx <- slope
   }
 
-  result_dims <- c(lon = dim(pr)[1] * nf,
-                   lat = dim(pr)[2] * nf,
-                   dim(pr)[3],
-                   member = nens)
-  names(result_dims)[3] <- names(dim(pr))[3]
+  result_dims <- c(dim(pr)[1] * nf, dim(pr)[2] * nf, dim(pr)[3], member = nens)
   r <- array(dim = result_dims)
 
   for (i in 1:nens) {
-- 
GitLab


From 57f8add4f14f2b5829c7b261cd3cd45f3644347c Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Tue, 19 Mar 2019 00:27:51 +0100
Subject: [PATCH 26/48] renamed stochastic ensemble dim to realization

---
 R/CST_RainFARM.R | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index d45c6ff4..58792bbd 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -158,7 +158,6 @@ RainFARM <- function(data, lon, lat, nf, weights = 1.,
   lat_f <- r[[2]]
 
   # Repeatedly apply .RainFARM
-  print(dim(data))
   result <- Apply(data, c("lon", "lat", time_dim), .RainFARM,
                   weights, nf, nens, slope, kmin,
                   fglob, fsmooth, verbose)$output1
@@ -180,7 +179,7 @@ RainFARM <- function(data, lon, lat, nf, weights = 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
-#' @return .RainFARM returns a downscaled array
+#' @return .RainFARM returns a downscaled array with dims (lon, lat, time, realization)
 #' @noRd
 
 .RainFARM <- function(pr, weights, nf, nens, slope, kmin,
@@ -194,7 +193,7 @@ RainFARM <- function(data, lon, lat, nf, weights = 1.,
     sx <- slope
   }
 
-  result_dims <- c(dim(pr)[1] * nf, dim(pr)[2] * nf, dim(pr)[3], member = nens)
+  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) {
-- 
GitLab


From 3f5b8be4edde9411d86efb96dbb7b379bccfed60 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Tue, 19 Mar 2019 01:54:50 +0100
Subject: [PATCH 27/48] reworked sorting of dims

---
 R/CST_RainFARM.R | 61 +++++++++++++++++++++++++++++++++---------------
 1 file changed, 42 insertions(+), 19 deletions(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index 58792bbd..eb6dead2 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -19,7 +19,7 @@
 #' @param fsmooth logical to conserve precipitation with a smnoothing kernel (default TRUE)
 #' @param time_dim name of time dimension (e.g. "ftime", "sdate", "time" ...)
 #' @param verbose logical for verbose output of the Julia package
-#' @param drop_ens_dim logical to remove the ensemble dimension if nens==1 (default FALSE)
+#' @param drop_realization_dim logical to remove the realization dimension if nens==1 (default TRUE)
 #'  
 #' @return CST_RainFARM() returns a downscaled CSTools object.
 #' @import multiApply
@@ -51,24 +51,44 @@ CST_RainFARM <- function(s2d_data, nf, weights = 1.,
                          slope = 0, kmin = 1, nens = 1,
                          fglob = FALSE, fsmooth = FALSE,
                          time_dim = NULL, verbose = FALSE,
-                         drop_ens_dim = FALSE) {
+                         drop_realization_dim = TRUE) {
 
+  # Check/detect time_dim (we need it also here)
+  if (is.null(time_dim)) {
+    time_dim_names <- c("ftime", "sdate", "time")
+    time_dim_num <- which(time_dim_names %in% names(dim(s2d_data$exp)))
+    if (length(time_dim_num) > 0) {
+      time_dim <- time_dim_names[time_dim_num[1]]
+    } else {
+      stop("Could not automatically detect a target time dimension.")
+    }
+  }
   res <- RainFARM(s2d_data$exp, s2d_data$lon, s2d_data$lat,
            nf, weights, nens,
            slope, kmin,
            fglob, fsmooth,
-           time_dim, drop_ens_dim,
+           time_dim, drop_realization_dim,
            verbose)
-  n_extra_dims <- length(dim(res$data)) - 3
-  extra_dims <- NULL
-  if (n_extra_dims > 0) {
-    extra_dims <- 1:n_extra_dims + 3
-  }
-  s2d_data$exp <- multiApply:::.aperm2(res$data, c(3, 2, 1, extra_dims))
 
+#  n_extra_dims <- length(dim(res$data)) - 3
+#  extra_dims <- NULL
+#  if (n_extra_dims > 0) {
+#    extra_dims <- 1:n_extra_dims + 3
+#  }
+
+# The input is already a s2d object, let's mantain the same order
+  if(length(dim(s2d_data$exp)) == length(dim(res$data))) {
+    iorder=sapply(names(dim(s2d_data$exp)), grep, names(dim(res$data)))
+  } else {
+    iorder=sapply(c(names(dim(s2d_data$exp)), "realization"), grep,
+                  names(dim(res$data)))
+  }
+  res$data=aperm(res$data, iorder)  
+# s2d_data$exp <- multiApply:::.aperm2(res$data, c(3, 2, 1, extra_dims))
+  s2d_data$exp <- res$data
   s2d_data$lon <- res$lon
   s2d_data$lat <- res$lat
-  s2d_data$obs <- NULL
+#  s2d_data$obs <- NULL 
 
   return(s2d_data)
 }
@@ -91,9 +111,9 @@ CST_RainFARM <- function(s2d_data, nf, weights = 1.,
 #' @param nens number of ensemble members to produce
 #' @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 time_dim name of time dimension (e.g. "ftime", "sdate", "time" ...)
+#' @param time_dim name of time dimension (e.g. "ftime", "sdate", "time"). If NULL (default)  one of these names is searched for. 
 #' @param verbose logical for verbose output of the Julia package
-#' @param drop_ens_dim logical to remove the ensemble dimension if nens==1
+#' @param drop_realization_dim logical to remove the ensemble dimension if nens==1
 #' @return RainFARM() returns a list containing the fine-scale longitudes, latitudes and the sequence of \code{nens} downscaled fields.
 #' @import JuliaCall
 #' @import multiApply
@@ -114,20 +134,20 @@ CST_RainFARM <- function(s2d_data, nf, weights = 1.,
 #' #     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)
 #' # downscale using weights (ww=1. means do not use weights)
-#' r <- RainFARM(pr, lon_mat, lat_mat, nf, fsmooth = TRUE, fglob = FALSE, weights = ww, nens = 2, verbose = TRUE)
+#' 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)
-#' #member  ftime	lat	lon
-#' # 	3 	20 	64 	64
+#' #  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, drop_ens_dim = FALSE,
+                     time_dim = NULL, drop_realization_dim = TRUE,
                      verbose = FALSE) {
 
   # Check/detect time_dim
@@ -162,15 +182,18 @@ RainFARM <- function(data, lon, lat, nf, weights = 1.,
                   weights, nf, nens, slope, kmin,
                   fglob, fsmooth, verbose)$output1
 
-  if (drop_ens_dim && nens == 1) {
-    dim(result) <- dim(result)[-which(names(dim(result)) == "member")[1]]
+  iorder=sapply(c(names(dim(data)), "realization"), grep, names(dim(result)))
+  result=aperm(result, iorder)  
+
+  if (drop_realization_dim && nens == 1) {
+    dim(result) <- dim(result)[-which(names(dim(result)) == "realization")[1]]
   }
 
   return(list(data = result, lon = lon_f, lat = lat_f))
 }
 
 #' Atomic RainFARM
-#' @param pr precipitation array to downscale
+#' @param pr precipitation array to downscale with dims (lon, lat, time)
 #' @param weights matrix (can be obtained using the \code{RFWeights} function (default 1. == none)
 #' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
 #' @param slope prescribed or automatic spectral slope (default 0. == automatic)
-- 
GitLab


From a894bae576c7cd7fc1e90f6a13909e257149f67c Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Tue, 19 Mar 2019 09:58:56 +0100
Subject: [PATCH 28/48] specify lon and lat names in RainFARM + better examples
 and docs

---
 R/CST_RainFARM.R | 59 ++++++++++++++++++------------------------------
 1 file changed, 22 insertions(+), 37 deletions(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index eb6dead2..02f37d74 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -36,39 +36,26 @@
 #' s2dv_data <- list(exp = exp, lon = lon, lat = lat)
 #' # Create a test array of weights
 #' ww <- array(1., dim = c(8 * nf, 8 * nf))
-#' res_s2dv <- CST_RainFARM(s2dv_data, nf, ww)
+#' res_s2dv <- CST_RainFARM(s2dv_data, nf, ww, nens=3)
 #' str(res_s2dv)
 #' #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_s2dv$exp)
-#' #  ftime     lat     lon dataset  member   sdate 
-#' #      4      64      64       1       2       3
+#' dataset  member   sdate   ftime     lat     lon    realization
+#'      1       2       3       4      64      64     3
 #'
 
-CST_RainFARM <- function(s2d_data, nf, weights = 1.,
-                         slope = 0, kmin = 1, nens = 1,
-                         fglob = FALSE, fsmooth = FALSE,
+CST_RainFARM <- function(s2d_data, nf, weights = 1., slope = 0, kmin = 1,
+                         nens = 1, fglob = FALSE, fsmooth = FALSE,
                          time_dim = NULL, verbose = FALSE,
                          drop_realization_dim = TRUE) {
 
-  # Check/detect time_dim (we need it also here)
-  if (is.null(time_dim)) {
-    time_dim_names <- c("ftime", "sdate", "time")
-    time_dim_num <- which(time_dim_names %in% names(dim(s2d_data$exp)))
-    if (length(time_dim_num) > 0) {
-      time_dim <- time_dim_names[time_dim_num[1]]
-    } else {
-      stop("Could not automatically detect a target time dimension.")
-    }
-  }
   res <- RainFARM(s2d_data$exp, s2d_data$lon, s2d_data$lat,
-           nf, weights, nens,
-           slope, kmin,
-           fglob, fsmooth,
-           time_dim, drop_realization_dim,
-           verbose)
+                  nf, weights, nens, slope, kmin, fglob, fsmooth,
+                  time_dim, lon_dim = "lon", lat_dim = "lat",
+                  drop_realization_dim, verbose)
 
 #  n_extra_dims <- length(dim(res$data)) - 3
 #  extra_dims <- NULL
@@ -76,19 +63,15 @@ CST_RainFARM <- function(s2d_data, nf, weights = 1.,
 #    extra_dims <- 1:n_extra_dims + 3
 #  }
 
-# The input is already a s2d object, let's mantain the same order
-  if(length(dim(s2d_data$exp)) == length(dim(res$data))) {
-    iorder=sapply(names(dim(s2d_data$exp)), grep, names(dim(res$data)))
-  } else {
-    iorder=sapply(c(names(dim(s2d_data$exp)), "realization"), grep,
-                  names(dim(res$data)))
-  }
-  res$data=aperm(res$data, iorder)  
+# The input is already a s2d object, let's mantain the same order 
+# RainFARM already returns an object with dims in the same order as input
+
 # s2d_data$exp <- multiApply:::.aperm2(res$data, c(3, 2, 1, extra_dims))
   s2d_data$exp <- res$data
   s2d_data$lon <- res$lon
   s2d_data$lat <- res$lat
-#  s2d_data$obs <- NULL 
+# Mantain what was in s2d_data
+# s2d_data$obs <- NULL 
 
   return(s2d_data)
 }
@@ -96,7 +79,8 @@ CST_RainFARM <- function(s2d_data, nf, weights = 1.,
 #' @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 an array with dims (lon, lat, time) in input.
+#' @description This function implements the RainFARM stochastic precipitation downscaling method
+#' and accepts an array with named dims ("lon", "lat") and a time dimension ("ftime", "sdate" or "time") in input.
 #' 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 julia_setup from the JuliaCall library, 
@@ -112,6 +96,8 @@ CST_RainFARM <- function(s2d_data, nf, weights = 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 time_dim name of time dimension (e.g. "ftime", "sdate", "time"). If NULL (default)  one of these names is searched for. 
+#' @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
 #' @param drop_realization_dim logical to remove the ensemble dimension if nens==1
 #' @return RainFARM() returns a list containing the fine-scale longitudes, latitudes and the sequence of \code{nens} downscaled fields.
@@ -144,11 +130,10 @@ CST_RainFARM <- function(s2d_data, nf, weights = 1.,
 #' #  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, drop_realization_dim = TRUE,
-                     verbose = FALSE) {
+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 = TRUE, verbose = FALSE) {
 
   # Check/detect time_dim
   if (is.null(time_dim)) {
@@ -178,7 +163,7 @@ RainFARM <- function(data, lon, lat, nf, weights = 1.,
   lat_f <- r[[2]]
 
   # Repeatedly apply .RainFARM
-  result <- Apply(data, c("lon", "lat", time_dim), .RainFARM,
+  result <- Apply(data, c(lon_dim, lat_dim, time_dim), .RainFARM,
                   weights, nf, nens, slope, kmin,
                   fglob, fsmooth, verbose)$output1
 
-- 
GitLab


From dabcf60a82da9f93fd04ce73b881f03786b1c108 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Tue, 19 Mar 2019 10:59:14 +0100
Subject: [PATCH 29/48] select a time dim only if length > 1

---
 R/CST_RainFARM.R | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index 02f37d74..bb7a5802 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -140,11 +140,18 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1,
     time_dim_names <- c("ftime", "sdate", "time")
     time_dim_num <- which(time_dim_names %in% names(dim(data)))
     if (length(time_dim_num) > 0) {
-      time_dim <- time_dim_names[time_dim_num[1]]
+    # Find time dimension with length > 1
+      ilong <- which(dim(data)[time_dim_names[time_dim_num]] > 1)
+      if (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
-- 
GitLab


From 7bc483840cf8097af73a911e2da16ea3b04157b1 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Tue, 19 Mar 2019 19:41:41 +0100
Subject: [PATCH 30/48] R/CST_RainFARM.R

---
 R/CST_RainFARM.R | 53 +++++++++++++++++++++++++++++-------------------
 1 file changed, 32 insertions(+), 21 deletions(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index bb7a5802..6b58f387 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -5,7 +5,8 @@
 #'
 #' @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.
-#' 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.
+#' 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 julia_setup from the JuliaCall library, 
 #' which takes a long time only the first time it is called on a machine.
 #'
@@ -13,11 +14,14 @@
 #' @param weights matrix (can be obtained using the \code{RFWeights} function (default 1. == none)
 #' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
 #' @param slope prescribed or automatic spectral slope (default 0. == automatic)
-#' @param kmin first wavenumber for automatic spectral slope
+#' @param kmin first wavenumber for spectral slope
 #' @param nens number of ensemble members to produce
 #' @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 time_dim name of time dimension (e.g. "ftime", "sdate", "time" ...)
+#' @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
 #' @param drop_realization_dim logical to remove the realization dimension if nens==1 (default TRUE)
 #'  
@@ -57,21 +61,9 @@ CST_RainFARM <- function(s2d_data, nf, weights = 1., slope = 0, kmin = 1,
                   time_dim, lon_dim = "lon", lat_dim = "lat",
                   drop_realization_dim, verbose)
 
-#  n_extra_dims <- length(dim(res$data)) - 3
-#  extra_dims <- NULL
-#  if (n_extra_dims > 0) {
-#    extra_dims <- 1:n_extra_dims + 3
-#  }
-
-# The input is already a s2d object, let's mantain the same order 
-# RainFARM already returns an object with dims in the same order as input
-
-# s2d_data$exp <- multiApply:::.aperm2(res$data, c(3, 2, 1, extra_dims))
   s2d_data$exp <- res$data
   s2d_data$lon <- res$lon
   s2d_data$lat <- res$lat
-# Mantain what was in s2d_data
-# s2d_data$obs <- NULL 
 
   return(s2d_data)
 }
@@ -82,7 +74,8 @@ CST_RainFARM <- function(s2d_data, nf, weights = 1., slope = 0, kmin = 1,
 #' @description This function implements the RainFARM stochastic precipitation downscaling method
 #' and accepts an array with named dims ("lon", "lat") and a time dimension ("ftime", "sdate" or "time") in input.
 #' 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.
+#' 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 julia_setup 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
@@ -91,11 +84,14 @@ CST_RainFARM <- function(s2d_data, nf, weights = 1., slope = 0, kmin = 1,
 #' @param weights matrix (can be obtained using the \code{RFWeights} function (default 1. == none)
 #' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
 #' @param slope prescribed or automatic spectral slope (default 0. == automatic)
-#' @param kmin first wavenumber for automatic spectral slope
+#' @param kmin first wavenumber for spectral slope
 #' @param nens number of ensemble members to produce
 #' @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 time_dim name of time dimension (e.g. "ftime", "sdate", "time"). If NULL (default)  one of these names is searched for. 
+#' @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
@@ -142,7 +138,7 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1,
     if (length(time_dim_num) > 0) {
     # Find time dimension with length > 1
       ilong <- which(dim(data)[time_dim_names[time_dim_num]] > 1)
-      if (ilong > 0) {
+      if (length(ilong) > 0) {
         time_dim <- time_dim_names[time_dim_num[ilong[1]]]
       } else {
         stop("No time dimension longer than one found.")
@@ -169,14 +165,29 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1,
   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, time_dim), .RainFARM,
+  result <- Apply(data, c(lon_dim, lat_dim, "rainfarm_samples"), .RainFARM,
                   weights, nf, nens, slope, kmin,
                   fglob, fsmooth, verbose)$output1
 
+  # Reorder as it was in data before apply + realization dim at end
   iorder=sapply(c(names(dim(data)), "realization"), grep, names(dim(result)))
   result=aperm(result, iorder)  
-
+  # Expand back to original dims + realization dim at end
+  dim(result) <- c(dim(result)[ind], cdim[-ind],
+                   dim(result)[length(dim(result))])
+  # Reorder as it was in data + realization dim at end
+  iorder=sapply(c(names(cdim0), "realization"), grep, names(dim(result)))
+  result=aperm(result, iorder)  
   if (drop_realization_dim && nens == 1) {
     dim(result) <- dim(result)[-which(names(dim(result)) == "realization")[1]]
   }
-- 
GitLab


From c3aae6fb27ef55733517abbe9653239426f8d844 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Tue, 19 Mar 2019 21:27:07 +0100
Subject: [PATCH 31/48] lint cleanup

---
 R/CST_RainFARM.R | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index 6b58f387..609a9978 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -168,26 +168,26 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1,
   # 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 <- 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]))  
- 
+  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
 
   # Reorder as it was in data before apply + realization dim at end
-  iorder=sapply(c(names(dim(data)), "realization"), grep, names(dim(result)))
-  result=aperm(result, iorder)  
+  iorder <- sapply(c(names(dim(data)), "realization"), grep, names(dim(result)))
+  result <- aperm(result, iorder)
   # Expand back to original dims + realization dim at end
   dim(result) <- c(dim(result)[ind], cdim[-ind],
                    dim(result)[length(dim(result))])
   # Reorder as it was in data + realization dim at end
-  iorder=sapply(c(names(cdim0), "realization"), grep, names(dim(result)))
-  result=aperm(result, iorder)  
+  iorder <- sapply(c(names(cdim0), "realization"), grep, names(dim(result)))
+  result <- aperm(result, iorder)
   if (drop_realization_dim && nens == 1) {
     dim(result) <- dim(result)[-which(names(dim(result)) == "realization")[1]]
   }
@@ -219,7 +219,8 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1,
     sx <- slope
   }
 
-  result_dims <- c(dim(pr)[1] * nf, dim(pr)[2] * nf, dim(pr)[3], realization = nens)
+  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) {
-- 
GitLab


From c1dd958e8866697539182d6410b21bffb986bc9e Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Tue, 19 Mar 2019 22:06:50 +0100
Subject: [PATCH 32/48] added RFSlope function to compute spectral slopes

---
 R/RFSlope.R | 146 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 146 insertions(+)
 create mode 100644 R/RFSlope.R

diff --git a/R/RFSlope.R b/R/RFSlope.R
new file mode 100644
index 00000000..41734c5f
--- /dev/null
+++ b/R/RFSlope.R
@@ -0,0 +1,146 @@
+#' @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 julia_setup from the JuliaCall library, 
+#' which takes a long time only the first time it is called on a machine.
+#'
+#' @param s2d_data CSTools object containing the precipitation array to downscale
+#' @param kmin first wavenumber for spectral slope (default 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 if A(k) is the spectral amplitude)
+#' @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))
+#' s2dv_data <- list(exp = exp, lon = lon, lat = lat)
+#' slopes <- CST_RFSlope(s2dv_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(s2d_data, kmin=1, time_dim=NULL) {
+
+  slopes <- RFSlope(s2d_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 julia_setup from the JuliaCall library, 
+#' which takes a long time only the first time it is called on a machine.
+#'
+#' @param data Array containing precipitation array to downscale
+#' @param kmin first wavenumber for spectral slope (default 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 if A(k) is the spectral amplitude)
+#' @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)
+#' slopes <- RFSlope(pr, kmin=1)
+#' 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 if 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)
+}
-- 
GitLab


From b916dcb7d9f17ca627c0252355803f7a082edbcb Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Tue, 19 Mar 2019 23:07:14 +0100
Subject: [PATCH 33/48] renamed rfslope

---
 R/{RFSlope.R => CST_RFSlope.R} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename R/{RFSlope.R => CST_RFSlope.R} (100%)

diff --git a/R/RFSlope.R b/R/CST_RFSlope.R
similarity index 100%
rename from R/RFSlope.R
rename to R/CST_RFSlope.R
-- 
GitLab


From 696d9c2b7a675959b0342f572fbd942fbddd39e1 Mon Sep 17 00:00:00 2001
From: Nicolau Manubens <nicolau.manubens@ic3.cat>
Date: Thu, 21 Mar 2019 20:26:49 +0100
Subject: [PATCH 34/48] Small modifications in documentation/parameter names.

---
 R/CST_RFSlope.R  | 12 ++++++------
 R/CST_RainFARM.R | 28 ++++++++++++++--------------
 R/RFWeights.R    |  2 +-
 3 files changed, 21 insertions(+), 21 deletions(-)

diff --git a/R/CST_RFSlope.R b/R/CST_RFSlope.R
index 41734c5f..38de21f3 100644
--- a/R/CST_RFSlope.R
+++ b/R/CST_RFSlope.R
@@ -8,7 +8,7 @@
 #' Please notice: The function uses julia_setup from the JuliaCall library, 
 #' which takes a long time only the first time it is called on a machine.
 #'
-#' @param s2d_data CSTools object containing the precipitation array to downscale
+#' @param data CSTools object containing the precipitation array to downscale
 #' @param kmin first wavenumber for spectral slope (default 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
@@ -26,8 +26,8 @@
 #' dim(lon) <- c(lon = length(lon))
 #' lat <- seq(40, 43.5, 0.5)
 #' dim(lat) <- c(lat = length(lat))
-#' s2dv_data <- list(exp = exp, lon = lon, lat = lat)
-#' slopes <- CST_RFSlope(s2dv_data)
+#' data <- list(exp = exp, lon = lon, lat = lat)
+#' slopes <- CST_RFSlope(data)
 #' dim(slopes)
 #' dataset  member   sdate
 #'      1       2       3 
@@ -36,9 +36,9 @@
 #' [1,] 1.893503 1.893503 1.893503
 #' [2,] 1.893503 1.893503 1.893503
 
-CST_RFSlope <- function(s2d_data, kmin=1, time_dim=NULL) {
+CST_RFSlope <- function(data, kmin = 1, time_dim = NULL) {
 
-  slopes <- RFSlope(s2d_data$exp, kmin, time_dim,
+  slopes <- RFSlope(data$exp, kmin, time_dim,
                      lon_dim = "lon", lat_dim = "lat")
 
   return(slopes)
@@ -127,7 +127,7 @@ RFSlope <- function(data, kmin = 1, time_dim = NULL,
   result <- Apply(data, c(lon_dim, lat_dim, "rainfarm_samples"),
                   .RFSlope, kmin)$output1
 
-  return(slopes=result)
+  return(slopes = result)
 }
 
 #' Atomic RFSlope
diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index 609a9978..b1390ab7 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -10,7 +10,7 @@
 #' Please notice: The function uses julia_setup from the JuliaCall library, 
 #' which takes a long time only the first time it is called on a machine.
 #'
-#' @param s2d_data CSTools object containing the precipitation array to downscale
+#' @param data CSTools object containing the precipitation array to downscale
 #' @param weights matrix (can be obtained using the \code{RFWeights} function (default 1. == none)
 #' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
 #' @param slope prescribed or automatic spectral slope (default 0. == automatic)
@@ -37,35 +37,35 @@
 #' dim(lon) <- c(lon = length(lon))
 #' lat <- seq(40, 43.5, 0.5)
 #' dim(lat) <- c(lat = length(lat))
-#' s2dv_data <- list(exp = exp, lon = lon, lat = 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_s2dv <- CST_RainFARM(s2dv_data, nf, ww, nens=3)
-#' str(res_s2dv)
+#' 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_s2dv$exp)
-#' dataset  member   sdate   ftime     lat     lon    realization
-#'      1       2       3       4      64      64     3
+#' dim(res$exp)
+#' # dataset  member   sdate   ftime     lat     lon    realization
+#' #       1       2       3       4      64      64     3
 #'
 
-CST_RainFARM <- function(s2d_data, nf, weights = 1., slope = 0, kmin = 1,
+CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1,
                          nens = 1, fglob = FALSE, fsmooth = FALSE,
                          time_dim = NULL, verbose = FALSE,
                          drop_realization_dim = TRUE) {
 
-  res <- RainFARM(s2d_data$exp, s2d_data$lon, s2d_data$lat,
+  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)
 
-  s2d_data$exp <- res$data
-  s2d_data$lon <- res$lon
-  s2d_data$lat <- res$lat
+  data$exp <- res$data
+  data$lon <- res$lon
+  data$lat <- res$lat
 
-  return(s2d_data)
+  return(data)
 }
 
 #' @rdname RainFARM
@@ -136,7 +136,7 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1,
     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
+      # 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]]]
diff --git a/R/RFWeights.R b/R/RFWeights.R
index 1d2de5a6..cd5ccb8d 100644
--- a/R/RFWeights.R
+++ b/R/RFWeights.R
@@ -59,7 +59,7 @@ RFWeights <- function(orofile, nf, lon=NULL, lat=NULL, reffile="",
   ww <- julia_call("rfweights", orofile, reffile, as.integer(nf),
                    weightsfn = weightsfn, varname = varname, fsmooth = fsmooth)
 
-  attributes(dim(ww))$names<-c("lon","lat")
+  attributes(dim(ww))$names <- c("lon","lat")
   unlink(reffile)
   return(ww)
 }
-- 
GitLab


From 7abab9ed11d38d614e4b3fa46084dbe475cf8126 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Fri, 22 Mar 2019 22:53:38 +0100
Subject: [PATCH 35/48] Upgraded documentation

---
 R/CST_RFSlope.R  |  66 ++++++++++++++--------
 R/CST_RainFARM.R | 140 +++++++++++++++++++++++++++++------------------
 R/RFWeights.R    |  53 +++++++++++-------
 3 files changed, 160 insertions(+), 99 deletions(-)

diff --git a/R/CST_RFSlope.R b/R/CST_RFSlope.R
index 41734c5f..d4f82825 100644
--- a/R/CST_RFSlope.R
+++ b/R/CST_RFSlope.R
@@ -1,21 +1,26 @@
 #' @rdname CST_RFSlope
 #' @title RainFARM spectral slopes from a CSTools object
 #'
-#' @author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
+#' @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 julia_setup from the JuliaCall library, 
+#' 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 s2d_data CSTools object containing the precipitation array to downscale
-#' @param kmin first wavenumber for spectral slope (default kmin=1)
-#' @param time_dim string or character array with name(s) of dimension(s) (e.g. "ftime", "sdate", "member" ...) 
+#' @param s2d_data CSTools object containing the spatial precipitation fields to downscale.
+#' The s2d_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 if A(k) is the spectral amplitude)
+#' @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
 #'
@@ -45,32 +50,45 @@ CST_RFSlope <- function(s2d_data, kmin=1, time_dim=NULL) {
 }
 
 #' @rdname RFSlope
-#' @title RainFARM Spectral Slopes from an array (reduced version)
+#' @title RainFARM spectral slopes from an array (reduced version)
 #'
-#' @author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
+#' @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 julia_setup from the JuliaCall library, 
+#' 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 precipitation array to downscale
-#' @param kmin first wavenumber for spectral slope (default 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
+#' @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 if A(k) is the spectral amplitude)
+#' 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
+#'
+#' # 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)
-#' slopes <- RFSlope(pr, kmin=1)
+#'
+#' # 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
@@ -131,10 +149,10 @@ RFSlope <- function(data, kmin = 1, time_dim = NULL,
 }
 
 #' Atomic RFSlope
-#' @param pr precipitation array to downscale with dims (lon, lat, time)
-#' @param kmin first wavenumber for spectral slope (default kmin=1)
+#' @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 if A(k) is the spectral amplitude)
+#' (the logarithmic slope of k*|A(k)|^2 where A(k) is the spectral amplitude).
 #' @noRd
 
 .RFSlope <- function(pr, kmin) {
diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index 609a9978..49e2409b 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -1,31 +1,47 @@
 #' @rdname CST_RainFARM
-#' @title RainFARM Stochastic Precipitation Downscaling of a CSTools object
+#' @title RainFARM stochastic precipitation downscaling of a CSTools object
 #'
-#' @author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
+#' @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.
-#' 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 julia_setup from the JuliaCall library, 
+#' @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 s2d_data CSTools object containing the precipitation array to downscale
-#' @param weights matrix (can be obtained using the \code{RFWeights} function (default 1. == none)
-#' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
-#' @param slope prescribed or automatic spectral slope (default 0. == automatic)
-#' @param kmin first wavenumber for spectral slope
-#' @param nens number of ensemble members to produce
-#' @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 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
+#' @param s2d_data CSTools object containing the spatial precipitation fields to downscale.
+#' The s2d_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
-#' @param drop_realization_dim logical to remove the realization dimension if nens==1 (default TRUE)
+#' 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 dimension
+#' if \code{nens==1} (default: TRUE).
 #'  
 #' @return CST_RainFARM() returns a downscaled CSTools object.
+#' If \code{nens>1} an additional dimension named "realization" is added to the \code{exp} array.
+#' The ordering of the dimensions in the \code{exp} element of the input object is mantained.
 #' @import multiApply
 #' @export
 #' @examples
@@ -69,34 +85,48 @@ CST_RainFARM <- function(s2d_data, nf, weights = 1., slope = 0, kmin = 1,
 }
 
 #' @rdname RainFARM
-#' @title RainFARM Stochastic Precipitation Downscaling (reduced version)
-#' @author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
+#' @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 an array with named dims ("lon", "lat") and a time dimension ("ftime", "sdate" or "time") in input.
+#' 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 julia_setup from the JuliaCall library, 
+#' 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
-#' @param lon vector or array of longitudes
-#' @param lat vector or array of latitudes
-#' @param weights matrix (can be obtained using the \code{RFWeights} function (default 1. == none)
-#' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
-#' @param slope prescribed or automatic spectral slope (default 0. == automatic)
-#' @param kmin first wavenumber for spectral slope
-#' @param nens number of ensemble members to produce
-#' @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 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
+#' @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
-#' @param drop_realization_dim logical to remove the ensemble dimension if nens==1
-#' @return RainFARM() returns a list containing the fine-scale longitudes, latitudes and the sequence of \code{nens} downscaled fields.
+#' 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 ensemble dimension if \code{nens==1} (default: TRUE).
+#' @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. The ordering of the input array dimensions is preserved.
 #' @import JuliaCall
 #' @import multiApply
 #' @export
@@ -196,16 +226,18 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1,
 }
 
 #' Atomic RainFARM
-#' @param pr precipitation array to downscale with dims (lon, lat, time)
-#' @param weights matrix (can be obtained using the \code{RFWeights} function (default 1. == none)
-#' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
-#' @param slope prescribed or automatic spectral slope (default 0. == automatic)
-#' @param kmin first wavenumber for automatic spectral slope
-#' @param nens number of ensemble members to produce
-#' @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
-#' @return .RainFARM returns a downscaled array with dims (lon, lat, time, realization)
+#' @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,
diff --git a/R/RFWeights.R b/R/RFWeights.R
index 1d2de5a6..f3f5ddff 100644
--- a/R/RFWeights.R
+++ b/R/RFWeights.R
@@ -1,39 +1,50 @@
-#' RainFARM Stochastic Precipitation Downscaling Weights
+#' Compute climatological weights for RainFARM stochastic precipitation downscaling
 #'
-#' @author Jost von Hardenberg, ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
+#' @author Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
 #'
-#' @description Compute 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 RFWeights() function uses julia_setup from the JuliaCall library, 
+#' @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 
-#' @param nf refinement factor for downscaling (the output resolution is increased by this factor)
-#' @param lon vector or array of longitudes (can be left out if \code{reffile} is used)
-#' @param lat vector or array of latitudes (can be left out 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 used
-#' @param varname optional name of variable to be read from reffile
-#' @param weightsfn optional filename to save the weights in netcdf format
-#' @param fsmooth logical to use smooth conservation (default) or large-scale box-average conservation
-#' @return The matrix containing the weights with dimensions (lon, lat).
+#' @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 CST_RainFARM() or RainFARM()
-#' # Specify a weightsfn filename if you wish to save the weights
-#' # Using an external fine-scale climatology file.
+#' # Create weights to be used with the CST_RainFARM() or RainFARM() functions
+#' # using an external fine-scale climatology file.
 #' 
 #' # 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", weightsfn="", fsmooth = TRUE)
+#' ww <- RFWeights("./worldclim.nc", nf, reffile="./input.nc")
 #' @export
 
 RFWeights <- function(orofile, nf, lon=NULL, lat=NULL, reffile="",
-- 
GitLab


From 957da5332d1b671b42d00b08df59617069700dbc Mon Sep 17 00:00:00 2001
From: Nicolau Manubens <nicolau.manubens@ic3.cat>
Date: Mon, 25 Mar 2019 14:49:48 +0100
Subject: [PATCH 36/48] Small fix in DESCRIPTION.

---
 DESCRIPTION | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index da0c26fa..1e096839 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -13,7 +13,7 @@ Authors@R: c(
     person("Veronica", "Torralba", , "veronica.torralba@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"),
+    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
-- 
GitLab


From fe3eabe093d5ae42898ee057fc6c6c6e03786ea2 Mon Sep 17 00:00:00 2001
From: Nicolau Manubens <nicolau.manubens@ic3.cat>
Date: Mon, 25 Mar 2019 14:50:36 +0100
Subject: [PATCH 37/48] Small fix in DESCRIPTION.

---
 DESCRIPTION | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 1e096839..eef988ee 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -26,7 +26,7 @@ Imports:
     abind,
     stats,
     JuliaCall,
-    ncdf4
+    ncdf4,
     ggplot2,
     plyr,
     data.table,
-- 
GitLab


From 0c3bc646a4b3e39b313cb69e72074b819706c6ad Mon Sep 17 00:00:00 2001
From: Nicolau Manubens <nicolau.manubens@ic3.cat>
Date: Mon, 25 Mar 2019 14:52:56 +0100
Subject: [PATCH 38/48] Small fix in DESCRIPTION.

---
 DESCRIPTION | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index eef988ee..ee5577e9 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -5,7 +5,7 @@ Version: 0.0.1
 Authors@R: c(
     person("BSC-CNS", role = c("aut", "cph")),
     person("Louis-Philippe", "Caron", , "louis-philippe.caron@bsc.es", role = "aut"),
-    person("Jost", "von Hardenberg", "j.vonhardenberg@isac.cnr.it", role = c("aut", "cph")),
+    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 = "aut"),
-- 
GitLab


From 290f2e0d1d7e71ffeafc5323ea2f9d6cb5f8306d Mon Sep 17 00:00:00 2001
From: Nicolau Manubens <nicolau.manubens@ic3.cat>
Date: Mon, 25 Mar 2019 15:16:11 +0100
Subject: [PATCH 39/48] Updated/generated RainFARM documentation.

---
 DESCRIPTION            |   2 +-
 NAMESPACE              |   7 +++
 R/CST_RainFARM.R       |  14 +++---
 man/CST_RFSlope.Rd     |  55 ++++++++++++++++++++++
 man/CST_RainFARM.Rd    |  89 ++++++++++++++++++++++++++++++++++++
 man/PlotForecastPDF.Rd |   9 ++--
 man/RFSlope.Rd         |  65 ++++++++++++++++++++++++++
 man/RFWeights.Rd       |  67 +++++++++++++++++++++++++++
 man/RainFARM.Rd        | 101 +++++++++++++++++++++++++++++++++++++++++
 9 files changed, 396 insertions(+), 13 deletions(-)
 create mode 100644 man/CST_RFSlope.Rd
 create mode 100644 man/CST_RainFARM.Rd
 create mode 100644 man/RFSlope.Rd
 create mode 100644 man/RFWeights.Rd
 create mode 100644 man/RainFARM.Rd

diff --git a/DESCRIPTION b/DESCRIPTION
index ee5577e9..a579abe5 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -37,4 +37,4 @@ Suggests:
 License: Apache License 2.0 | file LICENSE
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 6.0.1.9000
+RoxygenNote: 6.1.1
diff --git a/NAMESPACE b/NAMESPACE
index 286f913e..dad92f85 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_RainFARM.R b/R/CST_RainFARM.R
index b9511d4c..5679e188 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -42,8 +42,6 @@
 #' @return CST_RainFARM() returns a downscaled CSTools object.
 #' If \code{nens>1} an additional dimension named "realization" is added to the \code{exp} array.
 #' The ordering of the dimensions in the \code{exp} element of the input object is mantained.
-#' @import multiApply
-#' @export
 #' @examples
 #' #Example using CST_RainFARM for a CSTools object
 #' nf <- 8   # Choose a downscaling by factor 8
@@ -66,7 +64,8 @@
 #' # dataset  member   sdate   ftime     lat     lon    realization
 #' #       1       2       3       4      64      64     3
 #'
-
+#' @import multiApply
+#' @export
 CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1,
                          nens = 1, fglob = FALSE, fsmooth = FALSE,
                          time_dim = NULL, verbose = FALSE,
@@ -127,9 +126,6 @@ CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1,
 #' @param drop_realization_dim Logical to remove the ensemble dimension if \code{nens==1} (default: TRUE).
 #' @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. The ordering of the input array dimensions is preserved.
-#' @import JuliaCall
-#' @import multiApply
-#' @export
 #' @examples
 #' # Example for the 'reduced' RainFARM function 
 #' nf <- 8   # Choose a downscaling by factor 8
@@ -155,7 +151,10 @@ CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1,
 #' dim(res$data)
 #' #  lon         lat       ftime realization 
 #' #   64          64          20           2 
-
+#'
+#' @import JuliaCall
+#' @import multiApply
+#' @export
 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",
@@ -239,7 +238,6 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1,
 #' @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) {
 
diff --git a/man/CST_RFSlope.Rd b/man/CST_RFSlope.Rd
new file mode 100644
index 00000000..c97a4e75
--- /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 00000000..f2d66efb
--- /dev/null
+++ b/man/CST_RainFARM.Rd
@@ -0,0 +1,89 @@
+% 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 = FALSE, time_dim = NULL, verbose = FALSE,
+  drop_realization_dim = TRUE)
+}
+\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 dimension
+if \code{nens==1} (default: TRUE).}
+}
+\value{
+CST_RainFARM() returns a downscaled CSTools object.
+If \code{nens>1} an additional dimension named "realization" is added to the \code{exp} array.
+The ordering of the dimensions in the \code{exp} element of the input object is mantained.
+}
+\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 bd273aaf..0af93c07 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 00000000..650d5632
--- /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 00000000..61ebce18
--- /dev/null
+++ b/man/RFWeights.Rd
@@ -0,0 +1,67 @@
+% 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.
+
+# 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 00000000..0db8f8fa
--- /dev/null
+++ b/man/RainFARM.Rd
@@ -0,0 +1,101 @@
+% 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 = TRUE,
+  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 ensemble dimension if \code{nens==1} (default: TRUE).}
+
+\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. The ordering of the input array dimensions is preserved.
+}
+\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
+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}
+}
-- 
GitLab


From 138242d1a8bc16001157951f782dae3a6eaf62e3 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Tue, 26 Mar 2019 10:57:58 +0100
Subject: [PATCH 40/48] better treatment of realization dim

---
 R/CST_RainFARM.R | 35 +++++++++++++++++++++++------------
 1 file changed, 23 insertions(+), 12 deletions(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index 49e2409b..f8332fb4 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -70,7 +70,7 @@
 CST_RainFARM <- function(s2d_data, nf, weights = 1., slope = 0, kmin = 1,
                          nens = 1, fglob = FALSE, fsmooth = FALSE,
                          time_dim = NULL, verbose = FALSE,
-                         drop_realization_dim = TRUE) {
+                         drop_realization_dim = FALSE) {
 
   res <- RainFARM(s2d_data$exp, s2d_data$lon, s2d_data$lat,
                   nf, weights, nens, slope, kmin, fglob, fsmooth,
@@ -159,7 +159,7 @@ CST_RainFARM <- function(s2d_data, nf, weights = 1., slope = 0, kmin = 1,
 RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1,
                      slope = 0, kmin = 1, fglob = FALSE, fsmooth = TRUE,
                      time_dim = NULL, lon_dim = "lon", lat_dim = "lat",
-                     drop_realization_dim = TRUE, verbose = FALSE) {
+                     drop_realization_dim = FALSE, verbose = FALSE) {
 
   # Check/detect time_dim
   if (is.null(time_dim)) {
@@ -208,18 +208,29 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1,
   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 data before apply + realization dim at end
-  iorder <- sapply(c(names(dim(data)), "realization"), grep, names(dim(result)))
+  # 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)
-  # Expand back to original dims + realization dim at end
-  dim(result) <- c(dim(result)[ind], cdim[-ind],
-                   dim(result)[length(dim(result))])
-  # Reorder as it was in data + realization dim at end
-  iorder <- sapply(c(names(cdim0), "realization"), grep, names(dim(result)))
-  result <- aperm(result, iorder)
-  if (drop_realization_dim && nens == 1) {
-    dim(result) <- dim(result)[-which(names(dim(result)) == "realization")[1]]
+
+  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))
-- 
GitLab


From 4f6b56126536eb2009a3cd46c9c1560509cbff16 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Tue, 26 Mar 2019 11:25:09 +0100
Subject: [PATCH 41/48] doc ument new drop_realization_dim behaviour

---
 R/CST_RainFARM.R | 34 +++++++++++++++++++++++++++-------
 1 file changed, 27 insertions(+), 7 deletions(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index f8332fb4..681415b3 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -36,12 +36,20 @@
 #' 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 dimension
-#' if \code{nens==1} (default: TRUE).
-#'  
+#' @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.
-#' The ordering of the dimensions in the \code{exp} element of the input object is mantained.
+#' 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 mantained.
 #' @import multiApply
 #' @export
 #' @examples
@@ -124,9 +132,21 @@ CST_RainFARM <- function(s2d_data, nf, weights = 1., slope = 0, kmin = 1,
 #' @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 ensemble dimension if \code{nens==1} (default: TRUE).
+#' @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. The ordering of the input array dimensions is preserved.
+#' 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 mantained.
 #' @import JuliaCall
 #' @import multiApply
 #' @export
-- 
GitLab


From b95f35d7c20c479f9bd4c6f95961801aea6707ce Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Tue, 26 Mar 2019 11:34:44 +0100
Subject: [PATCH 42/48] typo

---
 R/CST_RainFARM.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index 681415b3..e7daeb77 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -49,7 +49,7 @@
 #' @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 mantained.
+#' The ordering of the remaining dimensions in the \code{exp} element of the input object is maintained.
 #' @import multiApply
 #' @export
 #' @examples
@@ -146,7 +146,7 @@ CST_RainFARM <- function(s2d_data, nf, weights = 1., slope = 0, kmin = 1,
 #' 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 mantained.
+#' The ordering of the remaining dimensions in the \code{exp} element of the input object is maintained.
 #' @import JuliaCall
 #' @import multiApply
 #' @export
-- 
GitLab


From 393a9ba121b825446181d748bc15fb89956f1fe7 Mon Sep 17 00:00:00 2001
From: Jost von Hardenberg <j.vonhardenberg@isac.cnr.it>
Date: Wed, 27 Mar 2019 00:08:27 +0100
Subject: [PATCH 43/48] fsmooth had wrong default value

---
 R/CST_RainFARM.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index e7daeb77..e957368d 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -76,7 +76,7 @@
 #'
 
 CST_RainFARM <- function(s2d_data, nf, weights = 1., slope = 0, kmin = 1,
-                         nens = 1, fglob = FALSE, fsmooth = FALSE,
+                         nens = 1, fglob = FALSE, fsmooth = TRUE,
                          time_dim = NULL, verbose = FALSE,
                          drop_realization_dim = FALSE) {
 
-- 
GitLab


From 20d6e813ab65bfeea15be263d1112d8cb5545f0a Mon Sep 17 00:00:00 2001
From: Nicolau Manubens <nicolau.manubens@ic3.cat>
Date: Wed, 27 Mar 2019 12:05:26 +0100
Subject: [PATCH 44/48] Small fixes.

---
 R/CST_RFSlope.R    | 20 ++++++++++----------
 R/CST_RainFARM.R   | 13 ++++++++-----
 man/CST_RFSlope.Rd | 10 +++++-----
 man/RFSlope.Rd     | 10 +++++-----
 man/RainFARM.Rd    | 13 ++++++++-----
 5 files changed, 36 insertions(+), 30 deletions(-)

diff --git a/R/CST_RFSlope.R b/R/CST_RFSlope.R
index 1155baec..c1e10db4 100644
--- a/R/CST_RFSlope.R
+++ b/R/CST_RFSlope.R
@@ -34,12 +34,12 @@
 #' data <- list(exp = exp, lon = lon, lat = lat)
 #' slopes <- CST_RFSlope(data)
 #' dim(slopes)
-#' dataset  member   sdate
-#'      1       2       3 
+#' # dataset  member   sdate
+#' #       1       2       3 
 #' slopes
-#'          [,1]     [,2]     [,3]
-#' [1,] 1.893503 1.893503 1.893503
-#' [2,] 1.893503 1.893503 1.893503
+#' #         [,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) {
 
@@ -93,11 +93,11 @@ CST_RFSlope <- function(data, kmin = 1, time_dim = NULL) {
 #' #  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
+#' #         [,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") {
diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index 5679e188..8162dba4 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -132,17 +132,20 @@ CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1,
 #' 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)
+#' 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)
+#' 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
-#' ww <- RFWeights("./worldclim.nc", nf, lon=lon_mat, lat=lat_mat, weightsfn="", fsmooth = TRUE)
+#' 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)
+#' 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 ...
diff --git a/man/CST_RFSlope.Rd b/man/CST_RFSlope.Rd
index c97a4e75..a9313740 100644
--- a/man/CST_RFSlope.Rd
+++ b/man/CST_RFSlope.Rd
@@ -43,12 +43,12 @@ 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 
+# dataset  member   sdate
+#       1       2       3 
 slopes
-         [,1]     [,2]     [,3]
-[1,] 1.893503 1.893503 1.893503
-[2,] 1.893503 1.893503 1.893503
+#         [,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/RFSlope.Rd b/man/RFSlope.Rd
index 650d5632..1f5689aa 100644
--- a/man/RFSlope.Rd
+++ b/man/RFSlope.Rd
@@ -54,11 +54,11 @@ 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
+#         [,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/RainFARM.Rd b/man/RainFARM.Rd
index 0db8f8fa..1cd7cabb 100644
--- a/man/RainFARM.Rd
+++ b/man/RainFARM.Rd
@@ -75,17 +75,20 @@ 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)
+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)
+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
-ww <- RFWeights("./worldclim.nc", nf, lon=lon_mat, lat=lat_mat, weightsfn="", fsmooth = TRUE)
+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)
+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 ...
-- 
GitLab


From c00ca1acc07bfbfdec9c35d30ae1cfafbcc847c8 Mon Sep 17 00:00:00 2001
From: Nicolau Manubens <nicolau.manubens@ic3.cat>
Date: Wed, 27 Mar 2019 12:49:29 +0100
Subject: [PATCH 45/48] Updates in the documentation.

---
 R/CST_RainFARM.R    |  5 +++++
 R/RFWeights.R       |  5 +++--
 man/CST_RainFARM.Rd | 20 ++++++++++++++------
 man/RFWeights.Rd    |  3 ++-
 man/RainFARM.Rd     | 18 +++++++++++++++---
 5 files changed, 39 insertions(+), 12 deletions(-)

diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index d912ef51..144cd67f 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -146,6 +146,9 @@ 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
 #' # Example for the 'reduced' RainFARM function 
 #' nf <- 8   # Choose a downscaling by factor 8
@@ -160,8 +163,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
+#' \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, 
diff --git a/R/RFWeights.R b/R/RFWeights.R
index be0b36e1..66b592e1 100644
--- a/R/RFWeights.R
+++ b/R/RFWeights.R
@@ -32,10 +32,10 @@
 #' @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)
@@ -45,6 +45,7 @@
 #'
 #' # 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="",
diff --git a/man/CST_RainFARM.Rd b/man/CST_RainFARM.Rd
index f2d66efb..2120194e 100644
--- a/man/CST_RainFARM.Rd
+++ b/man/CST_RainFARM.Rd
@@ -5,8 +5,8 @@
 \title{RainFARM stochastic precipitation downscaling of a CSTools object}
 \usage{
 CST_RainFARM(data, nf, weights = 1, slope = 0, kmin = 1, nens = 1,
-  fglob = FALSE, fsmooth = FALSE, time_dim = NULL, verbose = FALSE,
-  drop_realization_dim = TRUE)
+  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.
@@ -42,13 +42,21 @@ 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 dimension
-if \code{nens==1} (default: TRUE).}
+\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.
-The ordering of the dimensions in the \code{exp} element of the input object is mantained.
+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
diff --git a/man/RFWeights.Rd b/man/RFWeights.Rd
index 61ebce18..1f90924f 100644
--- a/man/RFWeights.Rd
+++ b/man/RFWeights.Rd
@@ -48,10 +48,10 @@ 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)
@@ -62,6 +62,7 @@ ww <- RFWeights("./worldclim.nc", nf, lon=lon_mat, lat=lat_mat, weightsfn="", fs
 # 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
index 1cd7cabb..b02b4102 100644
--- a/man/RainFARM.Rd
+++ b/man/RainFARM.Rd
@@ -6,7 +6,7 @@
 \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 = TRUE,
+  lon_dim = "lon", lat_dim = "lat", drop_realization_dim = FALSE,
   verbose = FALSE)
 }
 \arguments{
@@ -49,13 +49,23 @@ is searched and the first one with more than one element is chosen.}
 
 \item{lat_dim}{Name of lat dimension ("lat" by default).}
 
-\item{drop_realization_dim}{Logical to remove the ensemble dimension if \code{nens==1} (default: TRUE).}
+\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. The ordering of the input array dimensions is preserved.
+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
@@ -83,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
+\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, 
-- 
GitLab


From 040d328aaf7ec787b6af8428db1f6bd81682ccdd Mon Sep 17 00:00:00 2001
From: Nicolau Manubens <nicolau.manubens@ic3.cat>
Date: Wed, 27 Mar 2019 12:51:56 +0100
Subject: [PATCH 46/48] Updated GitLab CI with Julia module and R3.4 with
 JuliaCall.

---
 .gitlab-ci.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 48c15320..e234d4da 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -3,7 +3,7 @@ stages:
 build:
   stage: build
   script:
-    - module load R
+    - module load R/3.4.3-foss-2015a-bare Julia
     - R CMD build --resave-data .
     - R CMD check --as-cran CSTools_*.tar.gz
     - R -e 'covr::package_coverage()'
-- 
GitLab


From d936e44e4ff2e9d7bc2b2b9fa1b59e8d737a2f10 Mon Sep 17 00:00:00 2001
From: Nicolau Manubens <nicolau.manubens@ic3.cat>
Date: Wed, 27 Mar 2019 13:06:35 +0100
Subject: [PATCH 47/48] Fixes in GitLab CI.

---
 .gitlab-ci.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index e234d4da..b052e430 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -4,6 +4,7 @@ build:
   stage: build
   script:
     - 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 -e 'covr::package_coverage()'
-- 
GitLab


From 1be7c5c2e6605745ab6d41be33d49fad56824612 Mon Sep 17 00:00:00 2001
From: Nicolau Manubens <nicolau.manubens@ic3.cat>
Date: Wed, 27 Mar 2019 15:00:23 +0100
Subject: [PATCH 48/48] Fix in GitLab CI.

---
 .gitlab-ci.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index b052e430..ee44142d 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -6,5 +6,5 @@ build:
     - 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()'
-- 
GitLab