Merge pull request #69 from SantanderMetGroup/devel

merge devel to master

NAMESPACE

@@ -53,8 +53,10 @@ importFrom(stats,na.exclude)
 importFrom(stats,na.omit)
 importFrom(stats,nls)
 importFrom(stats,pgamma)
+importFrom(stats,pnorm)
 importFrom(stats,predict)
 importFrom(stats,qgamma)
+importFrom(stats,qnorm)
 importFrom(stats,quantile)
 importFrom(stats,resid)
 importFrom(stats,rgamma)

NEWS

@@ -2,6 +2,7 @@
 
 See the [Releases section](https://github.com/SantanderMetGroup/downscaleR/releases) for older version changes
 
+
 ## v3.0.0 (14 Jun 2018)
 
 * New user interface for flexible definition of predictors (`prepareData`) and prediction data (`prepareNewData`). 
@@ -49,3 +50,10 @@ See the [Releases section](https://github.com/SantanderMetGroup/downscaleR/relea
 * downscale.predict renamed --> downscalePredict
 * Other minor changes and documentation updates
 
+## v3.1.1 (5 Feb 2020)
+
+* Bug fix in `downscaleCV'
+* Improve the internal code of `downscaleChunk'
+* Add qdm and dqm bias correction methods
+* Bug fix in `biasCorrection'
+

R/downscaleCV.R

@@ -22,12 +22,13 @@
 #' @param y The observations dataset. It should be an object as returned by \pkg{loadeR}.
 #' @param method A string value. Type of transer function. Currently implemented options are \code{"analogs"}, \code{"GLM"} and \code{"NN"}.
 #' @param sampling.strategy Specifies a sampling strategy to define the training and test subsets. Possible values are 
-#' \code{"kfold.chronological"} (the default), \code{"kfold.random"} and \code{"leave-one-year-out"}.
+#' \code{"kfold.chronological"} (the default), \code{"kfold.random"}, \code{"leave-one-year-out"} and NULL.
 #' The \code{sampling.strategy} choices are next described:
 #' \itemize{
 #'   \item \code{"kfold.random"} creates the number of folds indicated in the \code{folds} argument by randomly sampling the entries along the time dimension.
 #'   \item \code{"kfold.chronological"} is similar to \code{"kfold.random"}, but the sampling is performed in ascending order along the time dimension.
-#'   \item \code{"leave-one-year-out"}. This schema performs a leave-one-year-out cross validation. It is equivalent to introduce in the argument \code{folds} a list of all years one by one.
+#'   \item \code{"leave-one-year-out"}. This scheme performs a leave-one-year-out cross validation. It is equivalent to introduce in the argument \code{folds} a list of all years one by one.
+#'   \item \code{NULL}. The folds are specified by the user in the function parameter \code{folds}.
 #' }
 #' The first two choices will be controlled by the argument \code{folds} (see below)
 #' @param folds This arguments controls the number of folds, or how these folds are created (ignored if \code{sampling.strategy = "leave-one-year-out"}). If it is given as a fraction in the range (0-1), 
@@ -44,7 +45,7 @@
 #' @param ... Optional parameters. These parameters are different depending on the method selected. 
 #' Every parameter has a default value set in the atomic functions in case that no selection is wanted. 
 #' Everything concerning these parameters is explained in the section \code{Details} of the function \code{\link[downscaleR]{downscaleTrain}}. However, if wanted, the atomic functions can be seen here: 
-#' \code{\link[downscaleR]{glm.train}} and \code{\link[deepnet]{nn.train}}.  
+#' \code{\link[downscaleR]{analogs.train}}, \code{\link[downscaleR]{glm.train}} and \code{\link[deepnet]{nn.train}}.  
 #' @details The function relies on \code{\link[downscaleR]{prepareData}}, \code{\link[downscaleR]{prepareNewData}}, \code{\link[downscaleR]{downscaleTrain}}, and \code{\link[downscaleR]{downscalePredict}}. 
 #' For more information please visit these functions. It is envisaged to allow for a flexible fine-tuning of the cross-validation scheme. It uses internally the \pkg{transformeR} 
 #' helper \code{\link[transformeR]{dataSplit}} for flexible data folding. 
@@ -120,24 +121,25 @@ downscaleCV <- function(x, y, method,
   x <- getTemporalIntersection(x,y,which.return = "obs")
   y <- getTemporalIntersection(x,y,which.return = "prd")
 
-  if (sampling.strategy == "leave-one-year-out") {
-    type <- "chronological"
-    folds <- as.list(getYearsAsINDEX(y) %>% unique())
-  }
-
-  if (sampling.strategy == "kfold.chronological") {
-    type <- "chronological"
-    if (!is.numeric(folds)) {
-      folds.user <- unlist(folds) %>% unique()
-      folds.data <- getYearsAsINDEX(y) %>% unique()
-      if (any(folds.user != folds.data)) stop("In the parameters folds you have indicated years that do not belong to the dataset. Please revise the setup of this parameter.")
+  if (!is.null(sampling.strategy)) {
+    if (sampling.strategy == "leave-one-year-out") {
+      type <- "chronological"
+      folds <- as.list(getYearsAsINDEX(y) %>% unique())
+    }
+
+    if (sampling.strategy == "kfold.chronological") {
+      type <- "chronological"
+      if (!is.numeric(folds)) {
+        folds.user <- unlist(folds) %>% unique() %>% sort()
+        folds.data <- getYearsAsINDEX(y) %>% unique()
+        if (any(folds.user != folds.data)) stop("In the parameters folds you have indicated years that do not belong to the dataset. Please revise the setup of this parameter.")
+      }
+    }
+    if (sampling.strategy == "kfold.random") {
+      type <- "random"
+      if (!is.numeric(folds)) stop("In kfold.random, the parameter folds represent the NUMBER of folds and thus, it should be a NUMERIC value.")
     }
   }
-  if (sampling.strategy == "kfold.random") {
-    type <- "random"
-    if (!is.numeric(folds)) stop("In kfold.random, the parameter folds represent the NUMBER of folds and thus, it should be a NUMERIC value.")
-  }
-
   if (is.list(folds)) {
     if (any(duplicated(unlist(folds)))) stop("Years can not appear in more than one fold")
   }

R/downscaleChunk.R

@@ -87,5 +87,13 @@ downscaleChunk <- function(x, y, newdata,
     })
     p <- NULL
   })
-  NULL
+
+  pred <- list()
+  for (i in 1:(length(newdata)+1)) {
+    lf <- list.files("./", pattern = paste0("dataset",i), full.names = TRUE)
+    chunk.list <- lapply(lf, function(x) get(load(x)))
+    pred[[i]] <- bindGrid(chunk.list, dimension = "lat")
+    file.remove(lf)
+  }
+  return(pred)
 }

R/isimip.R

@@ -79,9 +79,9 @@
 
 isimip <- function (y, x, newdata, threshold = 1, type = c("additive", "multiplicative")) {
 
-      obs <- y
-      pred <- x
-      sim <- newdata
+      obs <- redim(y, drop=TRUE)
+      pred <- redim(x, drop=TRUE)
+      sim <- redim(newdata, drop=TRUE)
 
       if (is.null(obs$Dates$start)){
             datesObs <- as.POSIXct(obs$Dates[[1]]$start, tz="GMT", format="%Y-%m-%d %H:%M:%S")
@@ -131,7 +131,6 @@ isimip <- function (y, x, newdata, threshold = 1, type = c("additive", "multipli
             callObs <- as.call(c(list(as.name("["),quote(pred$Data)), indTimeObs1))
             monthlyPred[indTimeObs] <- apply(eval(callObs), FUN = mean, MARGIN = setdiff(1:length(dimPred),pred.time.index), na.rm = TRUE)
       }
-
       if (is.null(sim$Dates$start)){
             datesFor <- as.POSIXct(sim$Dates[[1]]$start, tz="GMT", format="%Y-%m-%d %H:%M:%S")            
       }else{
@@ -294,8 +293,7 @@ isimip <- function (y, x, newdata, threshold = 1, type = c("additive", "multipli
                         }
                   }
             }
-      }
-      if(any(grepl(obs$Variable$varName,c("pr","tp","precipitation","precip")))){
+      } else if(any(grepl(obs$Variable$varName,c("pr","tp","precipitation","precip")))){
             if (length(threshold)==1){
                   threshold<-array(data = threshold, dim = 3)
             }
@@ -623,8 +621,7 @@ isimip <- function (y, x, newdata, threshold = 1, type = c("additive", "multipli
                   }
             }
             attr(sim$Data, "threshold") <- threshold
-      }
-      if(((any(grepl(obs$Variable$varName,c("radiation","pressure","wind","win dspeed","humidity","specific humidity","radiacion","presion","viento","humedad","humedad especifica","rss","rsds","rls","rlds","ps","slp","wss","huss","hus")))) | (type == "multiplicative")) & !multiField){
+      } else if(((any(grepl(obs$Variable$varName,c("radiation","pressure","wind","win dspeed","humidity","specific humidity","radiacion","presion","viento","humedad","humedad especifica","rss","rsds","rls","rlds","ps","slp","wss","huss","hus")))) | (type == "multiplicative")) & !multiField){
             if (length(threshold)==1){
                   threshold<-array(data = threshold, dim = 3)
             }
@@ -949,8 +946,7 @@ isimip <- function (y, x, newdata, threshold = 1, type = c("additive", "multipli
                   }
             }
             attr(sim$Data, "threshold") <- threshold
-      }
-      if((any(grepl(obs$Variable$varName,c("maximum temperature","temperatura maxima","tasmax","tmax","minimum temperature","temperatura minima","tasmin","tmin"))))){
+      } else if((any(grepl(obs$Variable$varName,c("maximum temperature","temperatura maxima","tasmax","tmax","minimum temperature","temperatura minima","tasmin","tmin"))))){
             indTas <- which(!is.na(match(obs$Variable$varName,c("tas","temperatura media","mean temperature","tmean"))))
             if (length(indTas) == 0){
                   stop("Mean temperature is needed to correct the Maximum and Minimum Temperatures")
@@ -1042,12 +1038,12 @@ isimip <- function (y, x, newdata, threshold = 1, type = c("additive", "multipli
                         }
                   }
             }
-      }
+            
       #       case {'uas';'vas';'ua';'va';'eastward wind component';'northward wind component'},
       #       if isempty(Ws),error('Wind speed is necessary for the correction of the eastward and northward wind component');end
       #       wsC=isimip(Ws.O,Ws.P,Ws.F,'variable','windspeed','datesobs',datesObs,'datesfor',datesFor);
       #       indC=find(~isnan(Ws.F) & Ws.F>0);F(indC)=(F(indC).*wsC(indC))./Ws.F(indC);
-      if((any(grepl(obs$Variable$varName,c("uas","vas","ua","va","eastward wind component","northward wind component"))))) {
+      } else if((any(grepl(obs$Variable$varName,c("uas","vas","ua","va","eastward wind component","northward wind component"))))) {
             indTas <- which(!is.na(match(obs$Variable$varName,c("wind","windspeed","viento","wss"))))
             if (length(indTas) == 0){
                   stop("Wind speed is needed to correct eastward and northward wind components")
@@ -1104,12 +1100,12 @@ isimip <- function (y, x, newdata, threshold = 1, type = c("additive", "multipli
                         sim$Data[indTasForAux[indCalmWind,]] <- calmWind
                   }
             }
-      }
+
       #       case {'prsn';'snowfall';'nieve'},
       #       if isempty(Pr),error('Precipitation is necessary for the correction of the snowfall');end
       #       prC=isimip(Pr.O,Pr.P,Pr.F,'variable','precipitation','datesobs',datesObs,'datesfor',datesFor,'threshold', threshold);
       #       indC=find(~isnan(Pr.F) & Pr.F>0);F(indC)=(F(indC).*prC(indC))./Pr.F(indC);
-      if((any(grepl(obs$Variable$varName,c("prsn","snowfall","nieve"))))) {
+      } else if((any(grepl(obs$Variable$varName,c("prsn","snowfall","nieve"))))) {
             indTas <- which(!is.na(match(obs$Variable$varName,c("pr","tp","precipitation","precip"))))
             if (length(indTas) == 0){
                   stop("Precipitation is needed to correct snow")

README.md

@@ -1,3 +1,5 @@
+[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3277316.svg)](https://doi.org/10.5281/zenodo.3277316)
+
 # What is downscaleR?
 
 **downscaleR** is an R package for empirical-statistical downscaling focusing on daily data and covering the most popular approaches (bias correction, Model Output Statistics, Perfect Prognosis) and techniques (e.g. quantile mapping, regression, analogs, neural networks). This package has been conceived to work in the framework of both seasonal forecasting and climate change studies. Thus, it considers ensemble members as a basic dimension of the data structure. Find out more about this package at the [downscaleR wiki](https://github.com/SantanderMetGroup/downscaleR/wiki).