From d3fff448b1ebae29d98b6bc9ded6af9c840f023a Mon Sep 17 00:00:00 2001 From: jbedia Date: Tue, 23 Jun 2015 10:33:16 +0200 Subject: [PATCH 1/5] Bug fix in aux date definition --- R/biasCorrection.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/biasCorrection.R b/R/biasCorrection.R index fdec185b..1e33fd8e 100644 --- a/R/biasCorrection.R +++ b/R/biasCorrection.R @@ -109,7 +109,7 @@ biasCorrection <- function (obs, pred, sim, method = c("qqmap", "delta", "scalin dimPerm[length(attr(obs$Data, "dimensions"))+k] <- indDiff } } - datesList <- as.POSIXct(obs$Dates$start, tz="GMT", format="%Y-%m-%d %H:%M:%S") + datesList <- as.POSIXct(obs$Dates$start, tz="GMT", format="%Y-%m-%d") yearList<-unlist(strsplit(as.character(datesList), "[-]")) dayListObs <- array(data = c(as.numeric(yearList[seq(2,length(yearList),3)]),as.numeric(yearList[seq(3,length(yearList),3)])), dim = c(length(datesList),2)) dayList <- unique(dayListObs,index.return = datesList) @@ -117,7 +117,7 @@ biasCorrection <- function (obs, pred, sim, method = c("qqmap", "delta", "scalin for (d in 1:dim(dayList)[1]){ indDays[which(sqrt((dayListObs[,1]-dayList[d,1])^2+(dayListObs[,2]-dayList[d,2])^2)==0)] <- d } - datesList <- as.POSIXct(sim$Dates$start, tz="GMT", format="%Y-%m-%d %H:%M:%S") + datesList <- as.POSIXct(sim$Dates$start, tz="GMT", format="%Y-%m-%d") yearList<-unlist(strsplit(as.character(datesList), "[-]")) dayListSim <- array(data = c(as.numeric(yearList[seq(2,length(yearList),3)]),as.numeric(yearList[seq(3,length(yearList),3)])), dim = c(length(datesList),2)) indDaysSim <- array(data = NaN, dim = c(length(datesList),1)) From cdd2eb313a12812702e41b077057c073729dbc32 Mon Sep 17 00:00:00 2001 From: miturbide Date: Fri, 26 Jun 2015 13:28:56 +0200 Subject: [PATCH 2/5] add gqm and qqmap extrapolation --- R/biasCorrection.R | 634 ++++++++++++++++++++++++++------------------- 1 file changed, 367 insertions(+), 267 deletions(-) diff --git a/R/biasCorrection.R b/R/biasCorrection.R index 1e33fd8e..5bfd344b 100644 --- a/R/biasCorrection.R +++ b/R/biasCorrection.R @@ -1,20 +1,25 @@ + #' @title Bias correction methods #' @description Implementation of several standard bias correction methods #' #' @template templateObsPredSim #' @param method method applied. Current accepted values are \code{"qqmap"}, \code{"delta"}, -#' \code{"scaling"}, \code{"unbiasing"} and \code{"piani"}. See details. +#' \code{"scaling"}, \code{"unbiasing"}, \code{"piani"} and \code{"gqm"}. See details. #' @param multi.member Should members be adjusted sepparately (TRUE, default), or jointly (FALSE)?. #' Ignored if the dataset has no members. #' @param pr.threshold The minimum value that is considered as a non-zero precipitation. Ignored for #' \code{varcode} values different from \code{"pr"}. Default to 1 (assuming mm). #' @param window Numeric value specifying the time window width used to calibrate. The window is centered on the target day. #' Default to \code{NULL}, which considers the whole period available. +#' @param extrapolation Character indicating the extrapolation method to be applied to correct values in +#' \code{"sim"} that are out of the range of \code{"pred"}. Extrapolation is applied only to the \code{"qqmap"} method, +#' thus, this argument is ignored if other bias correction method is selected. Default is \code{"no"} (do not extrapolate). +#' #' #' @details #' #' The methods available are \code{"qqmap"}, \code{"delta"}, \code{"unbiasing"}, -#' \code{"scaling"} and \code{"Piani"} (the latter used only for precipitation). +#' \code{"scaling"}, \code{"Piani"}, \code{"gqm"} (the two latter used only for precipitation). #' Next we make a brief description of each method: #' #' \strong{Delta} @@ -45,6 +50,12 @@ #' This method is described in Piani et al. 2010 and is applicable only to precipitation. It is based on the initial assumption that both observed #' and simulated intensity distributions are well approximated by the gamma distribution, therefore is a parametric q-q map #' that uses the theorical instead of the empirical distribution. +#' +#'\strong{Generalized Quantile Mapping (gqm)} +#' +#' This method is described in Gutjahr and Heinemann 2013. It is applicable only to precipitation and is similar to the Piani method. It applies a +#' gamma distribution to values under the threshold given by the 95th percentile (proosed by Yang et al. 2010) and a general Pareto +#' distribution (GPD) to values above the threshold. #' #' @seealso \code{\link{isimip}} for a trend-preserving method of model calibration #' @@ -60,29 +71,26 @@ #' \item A. Amengual, V. Homar, R. Romero, S. Alonso, and C. Ramis (2012) A Statistical Adjustment of Regional Climate Model Outputs to Local Scales: Application to Platja de Palma, Spain. J. Clim., 25, 939-957 #' #' \item C. Piani, J. O. Haerter and E. Coppola (2009) Statistical bias correction for daily precipitation in regional climate models over Europe, Theoretical and Applied Climatology, 99, 187-192 +#' +#' \item O. Gutjahr and G. Heinemann (2013) Comparing precipitation bias correction methods for high-resolution regional climate simulations using COSMO-CLM, Theoretical and Applied Climatology, 114, 511-529 #' } +#' +#' #' @author S. Herrera \email{sixto@@predictia.es} #' @export #' @examples \dontrun{ #' # These are the paths to the package built-in GSN and NCEP datasets -#' gsn.data.dir <- file.path(find.package("downscaleR.java"), -#' "datasets/observations/GSN_Iberia") -#' ncep.data.dir <- file.path(find.package("downscaleR.java"), -#' "datasets/reanalysis/Iberia_NCEP/Iberia_NCEP.ncml") +#' gsn.data.dir <- file.path(find.package("downscaleR"), "datasets/observations/GSN_Iberia") +#' ncep.data.dir <- file.path(find.package("downscaleR"), "datasets/reanalysis/Iberia_NCEP/Iberia_NCEP.ncml") #' gsn.inv <- dataInventory(gsn.data.dir) #' ncep.inv <- dataInventory(ncep.data.dir) #' str(gsn.inv) #' str(ncep.inv) -#' # Load precipitation for boreal winter (DJF) in the train (1991-2000) and test (2001-2010) periods, -#' # for the observations (GSN_Iberia) and the Iberia_NCEP datasets -#' obs <- loadStationData(dataset = gsn.data.dir, var="precip", lonLim = c(-12,10), latLim = c(33,47), -#' season=c(12,1,2), years = 1991:2000) -#' prd <- loadGridData(ncep.data.dir, var = "tp", lonLim = c(-12,10), latLim = c(33,47), -#' season = c(12,1,2), years = 1991:2000) -#' sim <- loadGridData(ncep.data.dir, var = "tp", lonLim = c(-12,10), latLim = c(33,47), -#' season = c(12,1,2), years = 2001:2010) -#' # Interpolation of the observations onto the grid of model: we use the method "nearest" -#' # and the getGrid function to ensure spatial consistency: +#' # Load precipitation for boreal winter (DJF) in the train (1991-2000) and test (2001-2010) periods, for the observations (GSN_Iberia) and the Iberia_NCEP datasets +#' obs <- loadStationData(dataset = gsn.data.dir, var="precip", lonLim = c(-12,10), latLim = c(33,47), season=c(12,1,2), years = 1991:2000) +#' prd <- loadGridData(ncep.data.dir, var = "tp", lonLim = c(-12,10), latLim = c(33,47), season = c(12,1,2), years = 1991:2000) +#' sim <- loadGridData(ncep.data.dir, var = "tp", lonLim = c(-12,10), latLim = c(33,47), season = c(12,1,2), years = 2001:2010) +#' # Interpolation of the observations onto the grid of model: we use the method "nearest" and the getGrid function to ensure spatial consistency: #' obs <- interpGridData(obs, new.grid = getGrid(prd), method = "nearest") #' # Apply the bias correction method: #' simBC <- biasCorrection (obs, prd, sim, method = "qqmap", pr.threshold = 1) # qq-mapping @@ -91,181 +99,182 @@ #' plotMeanField(simBC) #' } -biasCorrection <- function (obs, pred, sim, method = c("qqmap", "delta", "scaling", "unbiasing", "piani"), pr.threshold = 1, multi.member = TRUE, window = NULL) { - method <- match.arg(method, choices = c("qqmap", "delta", "unbiasing", "piani", "scaling")) - threshold <- pr.threshold - dimObs <- dim(obs$Data) - obs.time.index <- grep("^time$", attr(obs$Data, "dimensions")) - dimPred <- dim(pred$Data) - dimFor <- dim(sim$Data) - pred.time.index <- grep("^time$", attr(pred$Data, "dimensions")) - dimDiff <- NULL - dimPerm <- 1:length(attr(pred$Data, "dimensions")) - dimPerm[which(dimPerm<=length(dim(obs$Data)))] <- match(attr(obs$Data, "dimensions"), attr(pred$Data, "dimensions")) - if (length(setdiff(attr(pred$Data, "dimensions"),attr(obs$Data, "dimensions")))>0){ - dimDiff<-setdiff(attr(pred$Data, "dimensions"),attr(obs$Data, "dimensions")) - for (k in 1:length(dimDiff)) { - indDiff <- which(grepl(dimDiff[k],attr(pred$Data, "dimensions"))) - dimPerm[length(attr(obs$Data, "dimensions"))+k] <- indDiff - } - } - datesList <- as.POSIXct(obs$Dates$start, tz="GMT", format="%Y-%m-%d") - yearList<-unlist(strsplit(as.character(datesList), "[-]")) - dayListObs <- array(data = c(as.numeric(yearList[seq(2,length(yearList),3)]),as.numeric(yearList[seq(3,length(yearList),3)])), dim = c(length(datesList),2)) - dayList <- unique(dayListObs,index.return = datesList) - indDays <- array(data = NaN, dim = c(length(datesList),1)) - for (d in 1:dim(dayList)[1]){ - indDays[which(sqrt((dayListObs[,1]-dayList[d,1])^2+(dayListObs[,2]-dayList[d,2])^2)==0)] <- d - } - datesList <- as.POSIXct(sim$Dates$start, tz="GMT", format="%Y-%m-%d") - yearList<-unlist(strsplit(as.character(datesList), "[-]")) - dayListSim <- array(data = c(as.numeric(yearList[seq(2,length(yearList),3)]),as.numeric(yearList[seq(3,length(yearList),3)])), dim = c(length(datesList),2)) - indDaysSim <- array(data = NaN, dim = c(length(datesList),1)) - for (d in 1:dim(dayList)[1]){ - indDaysSim[which(sqrt((dayListSim[,1]-dayList[d,1])^2+(dayListSim[,2]-dayList[d,2])^2)==0)] <- d +biasCorrection <- function (obs, pred, sim, method = c("qqmap", "delta", "scaling", "unbiasing", "piani", "gqm"), pr.threshold = 1, multi.member = TRUE, window = NULL, extrapolation = c("no", "constant")) { + method <- match.arg(method, choices = c("qqmap", "delta", "unbiasing", "piani", "scaling", "gqm")) + threshold <- pr.threshold + dimObs <- dim(obs$Data) + obs.time.index <- grep("^time$", attr(obs$Data, "dimensions")) + dimPred <- dim(pred$Data) + dimFor <- dim(sim$Data) + pred.time.index <- grep("^time$", attr(pred$Data, "dimensions")) + dimDiff <- NULL + dimPerm <- 1:length(attr(pred$Data, "dimensions")) + dimPerm[which(dimPerm<=length(dim(obs$Data)))] <- match(attr(obs$Data, "dimensions"), attr(pred$Data, "dimensions")) + if (length(setdiff(attr(pred$Data, "dimensions"),attr(obs$Data, "dimensions")))>0){ + dimDiff<-setdiff(attr(pred$Data, "dimensions"),attr(obs$Data, "dimensions")) + for (k in 1:length(dimDiff)) { + indDiff <- which(grepl(dimDiff[k],attr(pred$Data, "dimensions"))) + dimPerm[length(attr(obs$Data, "dimensions"))+k] <- indDiff + } + } + datesList <- as.POSIXct(obs$Dates$start, tz="GMT", format="%Y-%m-%d") + yearList<-unlist(strsplit(as.character(datesList), "[-]")) + dayListObs <- array(data = c(as.numeric(yearList[seq(2,length(yearList),3)]),as.numeric(yearList[seq(3,length(yearList),3)])), dim = c(length(datesList),2)) + dayList <- unique(dayListObs,index.return = datesList) + indDays <- array(data = NaN, dim = c(length(datesList),1)) + for (d in 1:dim(dayList)[1]){ + indDays[which(sqrt((dayListObs[,1]-dayList[d,1])^2+(dayListObs[,2]-dayList[d,2])^2)==0)] <- d + } + datesList <- as.POSIXct(sim$Dates$start, tz="GMT", format="%Y-%m-%d") + yearList<-unlist(strsplit(as.character(datesList), "[-]")) + dayListSim <- array(data = c(as.numeric(yearList[seq(2,length(yearList),3)]),as.numeric(yearList[seq(3,length(yearList),3)])), dim = c(length(datesList),2)) + indDaysSim <- array(data = NaN, dim = c(length(datesList),1)) + for (d in 1:dim(dayList)[1]){ + indDaysSim[which(sqrt((dayListSim[,1]-dayList[d,1])^2+(dayListSim[,2]-dayList[d,2])^2)==0)] <- d + } + attrSim <- attr(sim$Data, "dimensions") + #apply function of calibration + if (length(dimDiff)==0){ + F <- calibrateProj(obs$Data, aperm(pred$Data,dimPerm), aperm(sim$Data,dimPerm), method = method, varcode = obs$Variable$varName, pr.threshold = threshold, extrapolate = extrapolation) + if (!any(dimPerm != 1:length(attr(pred$Data, "dimensions")))){ + sim$Data<-F + }else{ + sim$Data<-aperm(F,c(dimPerm[2:length(dim(obs$Data))],dimPerm[1])) + } + }else{ + if ((("member" %in% dimDiff) & isTRUE(multi.member)) | !("member" %in% dimDiff)) { + indDimDiff <- rep(list(bquote()), length(dimDiff)) + for (d in 1:length(dimDiff)){ + indDimDiff[[d]] <- 1:dimPred[match(dimDiff[d], attr(pred$Data, "dimensions"))] } - attrSim <- attr(sim$Data, "dimensions") - if (length(dimDiff)==0){ - F <- calibrateProj(obs$Data, aperm(pred$Data,dimPerm), aperm(sim$Data,dimPerm), method = method, varcode = obs$Variable$varName, pr.threshold = threshold) - if (!any(dimPerm != 1:length(attr(pred$Data, "dimensions")))){ - sim$Data<-F - }else{ - sim$Data<-aperm(F,c(dimPerm[2:length(dim(obs$Data))],dimPerm[1])) + indDimDiff<-as.matrix(expand.grid(indDimDiff)) + dimPermI <- match(attr(pred$Data, "dimensions")[setdiff(1:length(dimPred),match(dimDiff, attr(pred$Data, "dimensions")))],attr(obs$Data, "dimensions")) + for (i in 1:dim(indDimDiff)[1]){ + if (is.null(window)){ + indTimePrd <- rep(list(bquote()), length(dimPred)) + indTimeSim <- rep(list(bquote()), length(dimFor)) + indTimePrd[match(dimDiff, attr(pred$Data, "dimensions"))] <- as.list(indDimDiff[i,]) + indTimeSim[match(dimDiff, attr(sim$Data, "dimensions"))] <- as.list(indDimDiff[i,]) + callPrd <- as.call(c(list(as.name("["),quote(pred$Data)), indTimePrd)) + callSim <- as.call(c(list(as.name("["),quote(sim$Data)), indTimeSim)) + # attrSim <- attr(sim$Data, "dimensions") + F <- calibrateProj(aperm(obs$Data, dimPermI), eval(callPrd), eval(callSim), method = method, varcode = obs$Variable$varName, pr.threshold = threshold, extrapolate = extrapolation) + indTimeSim <- rep(list(bquote()), length(dimFor)) + for (d in 1:length(dimFor)){ + indTimeSim[[d]] <- 1:dimFor[d] + } + indTimeSim[match(dimDiff, attr(sim$Data, "dimensions"))] <- as.list(indDimDiff[i,]) + indTimeSim<-as.matrix(expand.grid(indTimeSim)) + sim$Data[indTimeSim] <- F + # attr(sim$Data, "dimensions") <- attrSim + }else{ + for (j in 1:dim(dayList)[1]){ + indObs <- which(indDays == j) + indObsWindow <- array(data = NaN, dim = c(window*length(indObs),1)) + for (d in 1:length(indObs)){ + indObsWindow[((d-1)*window+1):(d*window)] <- (indObs[d]-floor(window/2)):(indObs[d]+floor(window/2)) } - }else{ - if ((("member" %in% dimDiff) & isTRUE(multi.member)) | !("member" %in% dimDiff)) { - indDimDiff <- rep(list(bquote()), length(dimDiff)) - for (d in 1:length(dimDiff)){ - indDimDiff[[d]] <- 1:dimPred[match(dimDiff[d], attr(pred$Data, "dimensions"))] - } - indDimDiff<-as.matrix(expand.grid(indDimDiff)) - dimPermI <- match(attr(pred$Data, "dimensions")[setdiff(1:length(dimPred),match(dimDiff, attr(pred$Data, "dimensions")))],attr(obs$Data, "dimensions")) - for (i in 1:dim(indDimDiff)[1]){ - if (is.null(window)){ - indTimePrd <- rep(list(bquote()), length(dimPred)) - indTimeSim <- rep(list(bquote()), length(dimFor)) - indTimePrd[match(dimDiff, attr(pred$Data, "dimensions"))] <- as.list(indDimDiff[i,]) - indTimeSim[match(dimDiff, attr(sim$Data, "dimensions"))] <- as.list(indDimDiff[i,]) - callPrd <- as.call(c(list(as.name("["),quote(pred$Data)), indTimePrd)) - callSim <- as.call(c(list(as.name("["),quote(sim$Data)), indTimeSim)) -# attrSim <- attr(sim$Data, "dimensions") - F <- calibrateProj(aperm(obs$Data, dimPermI), eval(callPrd), eval(callSim), method = method, varcode = obs$Variable$varName, pr.threshold = threshold) - indTimeSim <- rep(list(bquote()), length(dimFor)) - for (d in 1:length(dimFor)){ - indTimeSim[[d]] <- 1:dimFor[d] - } - indTimeSim[match(dimDiff, attr(sim$Data, "dimensions"))] <- as.list(indDimDiff[i,]) - indTimeSim<-as.matrix(expand.grid(indTimeSim)) - sim$Data[indTimeSim] <- F -# attr(sim$Data, "dimensions") <- attrSim - }else{ - for (j in 1:dim(dayList)[1]){ - indObs <- which(indDays == j) - indObsWindow <- array(data = NaN, dim = c(window*length(indObs),1)) - for (d in 1:length(indObs)){ - indObsWindow[((d-1)*window+1):(d*window)] <- (indObs[d]-floor(window/2)):(indObs[d]+floor(window/2)) - } - indObsWindow[which(indObsWindow <= 0)] <- 1 - indObsWindow[which(indObsWindow > length(indDays))] <- length(indDays) - indObsWindow <- unique(indObsWindow) - indSim <- which(indDaysSim == j) - indTimeObs <- rep(list(bquote()), length(dimObs)) - indTimePrd <- rep(list(bquote()), length(dimPred)) - indTimeSim <- rep(list(bquote()), length(dimFor)) - # Consideramos solo la ventana temporal - indTimeObs[[obs.time.index]] <- indObsWindow - indTimePrd[[pred.time.index]] <- indObsWindow - indTimeSim[[pred.time.index]] <- indSim - indTimePrd[match(dimDiff, attr(pred$Data, "dimensions"))] <- as.list(indDimDiff[i,]) - indTimeSim[match(dimDiff, attr(sim$Data, "dimensions"))] <- as.list(indDimDiff[i,]) - callObs <- as.call(c(list(as.name("["),quote(obs$Data)), indTimeObs)) - callPrd <- as.call(c(list(as.name("["),quote(pred$Data)), indTimePrd)) - callSim <- as.call(c(list(as.name("["),quote(sim$Data)), indTimeSim)) - F <- calibrateProj(aperm(eval(callObs), dimPermI), eval(callPrd), eval(callSim), method = method, varcode = obs$Variable$varName, pr.threshold = threshold) - indTimeSim <- rep(list(bquote()), length(dimFor)) - for (d in 1:length(dimFor)){ - indTimeSim[[d]] <- 1:dimFor[d] - } - indTimeSim[[pred.time.index]] <- indSim - indTimeSim[match(dimDiff, attr(sim$Data, "dimensions"))] <- as.list(indDimDiff[i,]) - indTimeSim<-as.matrix(expand.grid(indTimeSim)) - sim$Data[indTimeSim] <- F - } -# attr(sim$Data, "dimensions") <- attrSim - } - } - }else{ - if (is.null(window)){ - obs.time.index <- grep("^time$", attr(obs$Data, "dimensions")) - pred.time.index <- grep("^time$", attr(pred$Data, "dimensions")) - pred.member.index <- grep("^member$", attr(pred$Data, "dimensions")) - auxPerm <- c(pred.time.index, pred.member.index,setdiff(1:length(dimPred),c(pred.time.index, pred.member.index))) - dimPermI <- match(attr(pred$Data, "dimensions")[setdiff(1:length(dimPred),match(dimDiff, attr(pred$Data, "dimensions")))],attr(obs$Data, "dimensions")) - auxPrd <- array(data = rep(aperm(pred$Data, auxPerm),1), dim = c(dimPred[pred.member.index]*dimPred[pred.time.index],dimPred[setdiff(1:length(dimPred),c(pred.time.index, pred.member.index))])) - auxSim <- array(data = rep(aperm(sim$Data, auxPerm),1), dim = c(dimFor[pred.member.index]*dimFor[pred.time.index],dimFor[setdiff(1:length(dimFor),c(pred.time.index, pred.member.index))])) - auxObs <- array(data = NaN, dim = c(dimPred[pred.member.index]*dimPred[pred.time.index],dimPred[setdiff(1:length(dimPred),c(pred.time.index, pred.member.index))])) - for (i in 1:dimPred[pred.member.index]){ - indMember <- ((i-1)*dimPred[pred.time.index]+1):(i*dimPred[pred.time.index]) - auxObs[indMember,,] <- aperm(obs$Data, dimPermI) - } - # attrSim <- attr(sim$Data, "dimensions") - F <- calibrateProj(auxObs, auxPrd, auxSim, method = method, varcode = obs$Variable$varName, pr.threshold = threshold) - sim$Data <- aperm(array(data = rep(F,1), dim = c(dimFor[pred.time.index],dimFor[pred.member.index],dimFor[setdiff(1:length(dimFor),c(pred.time.index, pred.member.index))])), auxPerm) -# attr(sim$Data, "dimensions") <- attrSim - }else{ - obs.time.index <- grep("^time$", attr(obs$Data, "dimensions")) - pred.time.index <- grep("^time$", attr(pred$Data, "dimensions")) - pred.member.index <- grep("^member$", attr(pred$Data, "dimensions")) - auxPerm <- c(pred.time.index, pred.member.index,setdiff(1:length(dimPred),c(pred.time.index, pred.member.index))) - dimPermI <- match(attr(pred$Data, "dimensions")[setdiff(1:length(dimPred),match(dimDiff, attr(pred$Data, "dimensions")))],attr(obs$Data, "dimensions")) - for (j in 1:dim(dayList)[1]){ - indObs <- which(indDays == j) - indObsWindow <- array(data = NaN, dim = c(window*length(indObs),1)) - for (d in 1:length(indObs)){ - indObsWindow[((d-1)*window+1):(d*window)] <- (indObs[d]-floor(window/2)):(indObs[d]+floor(window/2)) - } - indObsWindow[which(indObsWindow <= 0)] <- 1 - indObsWindow[which(indObsWindow > length(indDays))] <- length(indDays) - indObsWindow <- unique(indObsWindow) - indSim <- which(indDaysSim == j) - indTimeObs <- rep(list(bquote()), length(dimObs)) - indTimePrd <- rep(list(bquote()), length(dimPred)) - indTimeSim <- rep(list(bquote()), length(dimFor)) - # Consideramos solo la ventana temporal - indTimeObs[[obs.time.index]] <- indObsWindow - indTimePrd[[pred.time.index]] <- indObsWindow - indTimeSim[[pred.time.index]] <- indSim - callObs <- as.call(c(list(as.name("["),quote(obs$Data)), indTimeObs)) - callPrd <- as.call(c(list(as.name("["),quote(pred$Data)), indTimePrd)) - callSim <- as.call(c(list(as.name("["),quote(sim$Data)), indTimeSim)) - auxPrd <- array(data = rep(aperm(eval(callPrd), auxPerm),1), dim = c(dimPred[pred.member.index]*length(indObsWindow),dimPred[setdiff(1:length(dimPred),c(pred.time.index, pred.member.index))])) - auxSim <- array(data = rep(aperm(eval(callSim), auxPerm),1), dim = c(dimFor[pred.member.index]*length(indSim),dimFor[setdiff(1:length(dimFor),c(pred.time.index, pred.member.index))])) - auxObs <- array(data = NaN, dim = c(dimPred[pred.member.index]*length(indObsWindow),dimPred[setdiff(1:length(dimPred),c(pred.time.index, pred.member.index))])) - for (i in 1:dimPred[pred.member.index]){ - indMember <- ((i-1)*length(indObsWindow)+1):(i*length(indObsWindow)) - auxObs[indMember,,] <- aperm(eval(callObs), dimPermI) - } - # attrSim <- attr(sim$Data, "dimensions") - F <- calibrateProj(auxObs, auxPrd, auxSim, method = method, varcode = obs$Variable$varName, pr.threshold = threshold) - F <- aperm(array(data = rep(F,1), dim = c(length(indSim),dimFor[pred.member.index],dimFor[setdiff(1:length(dimFor),c(pred.time.index, pred.member.index))])), auxPerm) - indTimeSim <- rep(list(bquote()), length(dimFor)) - for (d in 1:length(dimFor)){ - indTimeSim[[d]] <- 1:dimFor[d] - } - indTimeSim[[pred.time.index]] <- indSim - indTimeSim<-as.matrix(expand.grid(indTimeSim)) - sim$Data[indTimeSim] <- F -# attr(sim$Data, "dimensions") <- attrSim - } - } + indObsWindow[which(indObsWindow <= 0)] <- 1 + indObsWindow[which(indObsWindow > length(indDays))] <- length(indDays) + indObsWindow <- unique(indObsWindow) + indSim <- which(indDaysSim == j) + indTimeObs <- rep(list(bquote()), length(dimObs)) + indTimePrd <- rep(list(bquote()), length(dimPred)) + indTimeSim <- rep(list(bquote()), length(dimFor)) + # Consideramos solo la ventana temporal + indTimeObs[[obs.time.index]] <- indObsWindow + indTimePrd[[pred.time.index]] <- indObsWindow + indTimeSim[[pred.time.index]] <- indSim + indTimePrd[match(dimDiff, attr(pred$Data, "dimensions"))] <- as.list(indDimDiff[i,]) + indTimeSim[match(dimDiff, attr(sim$Data, "dimensions"))] <- as.list(indDimDiff[i,]) + callObs <- as.call(c(list(as.name("["),quote(obs$Data)), indTimeObs)) + callPrd <- as.call(c(list(as.name("["),quote(pred$Data)), indTimePrd)) + callSim <- as.call(c(list(as.name("["),quote(sim$Data)), indTimeSim)) + F <- calibrateProj(aperm(eval(callObs), dimPermI), eval(callPrd), eval(callSim), method = method, varcode = obs$Variable$varName, pr.threshold = threshold, extrapolate = extrapolation) + indTimeSim <- rep(list(bquote()), length(dimFor)) + for (d in 1:length(dimFor)){ + indTimeSim[[d]] <- 1:dimFor[d] } + indTimeSim[[pred.time.index]] <- indSim + indTimeSim[match(dimDiff, attr(sim$Data, "dimensions"))] <- as.list(indDimDiff[i,]) + indTimeSim<-as.matrix(expand.grid(indTimeSim)) + sim$Data[indTimeSim] <- F + } + # attr(sim$Data, "dimensions") <- attrSim + } } - if (any(grepl(obs$Variable$varName,c("pr","tp","precipitation","precip")))){ - attr(sim$Data, "threshold") <- threshold + }else{ + if (is.null(window)){ + obs.time.index <- grep("^time$", attr(obs$Data, "dimensions")) + pred.time.index <- grep("^time$", attr(pred$Data, "dimensions")) + pred.member.index <- grep("^member$", attr(pred$Data, "dimensions")) + auxPerm <- c(pred.time.index, pred.member.index,setdiff(1:length(dimPred),c(pred.time.index, pred.member.index))) + dimPermI <- match(attr(pred$Data, "dimensions")[setdiff(1:length(dimPred),match(dimDiff, attr(pred$Data, "dimensions")))],attr(obs$Data, "dimensions")) + auxPrd <- array(data = rep(aperm(pred$Data, auxPerm),1), dim = c(dimPred[pred.member.index]*dimPred[pred.time.index],dimPred[setdiff(1:length(dimPred),c(pred.time.index, pred.member.index))])) + auxSim <- array(data = rep(aperm(sim$Data, auxPerm),1), dim = c(dimFor[pred.member.index]*dimFor[pred.time.index],dimFor[setdiff(1:length(dimFor),c(pred.time.index, pred.member.index))])) + auxObs <- array(data = NaN, dim = c(dimPred[pred.member.index]*dimPred[pred.time.index],dimPred[setdiff(1:length(dimPred),c(pred.time.index, pred.member.index))])) + for (i in 1:dimPred[pred.member.index]){ + indMember <- ((i-1)*dimPred[pred.time.index]+1):(i*dimPred[pred.time.index]) + auxObs[indMember,,] <- aperm(obs$Data, dimPermI) + } + # attrSim <- attr(sim$Data, "dimensions") + F <- calibrateProj(auxObs, auxPrd, auxSim, method = method, varcode = obs$Variable$varName, pr.threshold = threshold, extrapolate = extrapolation) + sim$Data <- aperm(array(data = rep(F,1), dim = c(dimFor[pred.time.index],dimFor[pred.member.index],dimFor[setdiff(1:length(dimFor),c(pred.time.index, pred.member.index))])), auxPerm) + # attr(sim$Data, "dimensions") <- attrSim + }else{ + obs.time.index <- grep("^time$", attr(obs$Data, "dimensions")) + pred.time.index <- grep("^time$", attr(pred$Data, "dimensions")) + pred.member.index <- grep("^member$", attr(pred$Data, "dimensions")) + auxPerm <- c(pred.time.index, pred.member.index,setdiff(1:length(dimPred),c(pred.time.index, pred.member.index))) + dimPermI <- match(attr(pred$Data, "dimensions")[setdiff(1:length(dimPred),match(dimDiff, attr(pred$Data, "dimensions")))],attr(obs$Data, "dimensions")) + for (j in 1:dim(dayList)[1]){ + indObs <- which(indDays == j) + indObsWindow <- array(data = NaN, dim = c(window*length(indObs),1)) + for (d in 1:length(indObs)){ + indObsWindow[((d-1)*window+1):(d*window)] <- (indObs[d]-floor(window/2)):(indObs[d]+floor(window/2)) + } + indObsWindow[which(indObsWindow <= 0)] <- 1 + indObsWindow[which(indObsWindow > length(indDays))] <- length(indDays) + indObsWindow <- unique(indObsWindow) + indSim <- which(indDaysSim == j) + indTimeObs <- rep(list(bquote()), length(dimObs)) + indTimePrd <- rep(list(bquote()), length(dimPred)) + indTimeSim <- rep(list(bquote()), length(dimFor)) + # Consideramos solo la ventana temporal + indTimeObs[[obs.time.index]] <- indObsWindow + indTimePrd[[pred.time.index]] <- indObsWindow + indTimeSim[[pred.time.index]] <- indSim + callObs <- as.call(c(list(as.name("["),quote(obs$Data)), indTimeObs)) + callPrd <- as.call(c(list(as.name("["),quote(pred$Data)), indTimePrd)) + callSim <- as.call(c(list(as.name("["),quote(sim$Data)), indTimeSim)) + auxPrd <- array(data = rep(aperm(eval(callPrd), auxPerm),1), dim = c(dimPred[pred.member.index]*length(indObsWindow),dimPred[setdiff(1:length(dimPred),c(pred.time.index, pred.member.index))])) + auxSim <- array(data = rep(aperm(eval(callSim), auxPerm),1), dim = c(dimFor[pred.member.index]*length(indSim),dimFor[setdiff(1:length(dimFor),c(pred.time.index, pred.member.index))])) + auxObs <- array(data = NaN, dim = c(dimPred[pred.member.index]*length(indObsWindow),dimPred[setdiff(1:length(dimPred),c(pred.time.index, pred.member.index))])) + for (i in 1:dimPred[pred.member.index]){ + indMember <- ((i-1)*length(indObsWindow)+1):(i*length(indObsWindow)) + auxObs[indMember,,] <- aperm(eval(callObs), dimPermI) + } + # attrSim <- attr(sim$Data, "dimensions") + F <- calibrateProj(auxObs, auxPrd, auxSim, method = method, varcode = obs$Variable$varName, pr.threshold = threshold, extrapolate = extrapolation) + F <- aperm(array(data = rep(F,1), dim = c(length(indSim),dimFor[pred.member.index],dimFor[setdiff(1:length(dimFor),c(pred.time.index, pred.member.index))])), auxPerm) + indTimeSim <- rep(list(bquote()), length(dimFor)) + for (d in 1:length(dimFor)){ + indTimeSim[[d]] <- 1:dimFor[d] + } + indTimeSim[[pred.time.index]] <- indSim + indTimeSim<-as.matrix(expand.grid(indTimeSim)) + sim$Data[indTimeSim] <- F + # attr(sim$Data, "dimensions") <- attrSim + } } - attr(sim$Data, "dimensions") <- attrSim - attr(sim$Data, "correction") <- method - return(sim) + } + } + if (any(grepl(obs$Variable$varName,c("pr","tp","precipitation","precip")))){ + attr(sim$Data, "threshold") <- threshold + } + attr(sim$Data, "dimensions") <- attrSim + attr(sim$Data, "correction") <- method + return(sim) } # End ################################################################################ @@ -276,7 +285,7 @@ biasCorrection <- function (obs, pred, sim, method = c("qqmap", "delta", "scalin #' #' @template templateObsPredSim #' @param method method applied. Current accepted values are \code{"qqmap"}, \code{"delta"}, -#' \code{"scaling"}, \code{"unbiasing"} and \code{"piani"}. See details. +#' \code{"scaling"}, \code{"unbiasing"}, \code{"piani"} and \code{"gqm"}. See details. #' @param varcode Variable code. This is not the variable itself to be corrected, but #' rather it referes to its nature and distributional properties. For instance, \code{"tas"} applies for #' temperature-like variables (i.e.: unbounded gaussian variables that can take both negative and positive values), @@ -285,14 +294,21 @@ biasCorrection <- function (obs, pred, sim, method = c("qqmap", "delta", "scalin #' but without an upper bound) and \code{"pr"}, specifically for precipitation. #' @param pr.threshold The minimum value that is considered as a non-zero precipitation. Ignored for #' \code{varcode} values different from \code{"pr"}. Default to 1 (assuming mm). -#' @author S. Herrera \email{sixto@@predictia.es} +#' @param extrapolate Character indicating the extrapolation method to be applied to correct values in +#' \code{"sim"} that are out of the range of \code{"pred"}. Extrapolation is applied only to the \code{"qqmap"} method, +#' thus, this argument is ignored if other bias correction method is selected. Default is \code{"no"} (do not extrapolate). +#' +#@author S. Herrera \email{sixto@@predictia.es} #' @export #' @family downscaling #' @family calibration #' @importFrom MASS fitdistr +#' @importFrom evd fpot #' @keywords internal #' -calibrateProj <- function (obs, pred, sim, method = c("qqmap", "delta", "scaling", "unbiasing", "piani"), varcode = c("tas", "hurs", "tp", "pr", "wss"), pr.threshold = 1) { + + +calibrateProj <- function (obs, pred, sim, method = c("qqmap", "delta", "scaling", "unbiasing", "piani", "gqm"), varcode = c("tas", "hurs", "tp", "pr", "wss"), pr.threshold = 1, extrapolate = c("no", "constant")) { if (any(grepl(varcode,c("pr","tp","precipitation","precip")))) { threshold<-pr.threshold nP<-matrix(data = NA, ncol=dim(pred)[3], nrow=dim(pred)[2]) @@ -308,10 +324,10 @@ calibrateProj <- function (obs, pred, sim, method = c("qqmap", "delta", "scaling Os<-sort(obs[,i,j], decreasing = FALSE, na.last = NA) ind<-which(Ps > threshold & !is.na(Ps)) if (length(ind)==0){ - ind <- max(which(!is.na(Ps))) - ind <- min(c(length(Os),ind)) + ind <- max(which(!is.na(Ps))) + ind <- min(c(length(Os),ind)) }else{ - ind<-min(which(Ps > threshold & !is.na(Ps))) + ind<-min(which(Ps > threshold & !is.na(Ps))) } # [Shape parameter Scale parameter] if (length(unique(Os[(nP[i,j]+1):ind]))<6){ @@ -339,88 +355,172 @@ calibrateProj <- function (obs, pred, sim, method = c("qqmap", "delta", "scaling } } if (method == "delta") { - if (dim(sim)[1]!=dim(obs)[1]){ - stop("sim and obs should have the same dimensions") - }else{ - simMean <- apply(sim, FUN = mean, MARGIN = c(2,3), na.rm = TRUE) - prdMean <- apply(pred, FUN = mean, MARGIN = c(2,3), na.rm = TRUE) - for (k in 1:min(c(dim(sim)[1],dim(obs)[1]))) { - sim[k,,]<-obs[k,,]-prdMean+simMean - } - } - } - if (method == "unbiasing") { - obsMean <- apply(obs, FUN = mean, MARGIN = c(2,3), na.rm = TRUE) - prdMean <- apply(pred, FUN = mean, MARGIN = c(2,3), na.rm = TRUE) - for (k in 1:dim(sim)[1]) { - sim[k,,]<-sim[k,,]-prdMean+obsMean - } + if (dim(sim)[1]!=dim(obs)[1]){ + stop("sim and obs should have the same dimensions") + }else{ + simMean <- apply(sim, FUN = mean, MARGIN = c(2,3), na.rm = TRUE) + prdMean <- apply(pred, FUN = mean, MARGIN = c(2,3), na.rm = TRUE) + for (k in 1:min(c(dim(sim)[1],dim(obs)[1]))) { + sim[k,,]<-obs[k,,]-prdMean+simMean } - if (method == "scaling") { - obsMean <- apply(obs, FUN = mean, MARGIN = c(2,3), na.rm = TRUE) - prdMean <- apply(pred, FUN = mean, MARGIN = c(2,3), na.rm = TRUE) - for (k in 1:dim(sim)[1]) { - sim[k,,]<-(sim[k,,]/prdMean)*obsMean + } + } + if (method == "unbiasing") { + obsMean <- apply(obs, FUN = mean, MARGIN = c(2,3), na.rm = TRUE) + prdMean <- apply(pred, FUN = mean, MARGIN = c(2,3), na.rm = TRUE) + for (k in 1:dim(sim)[1]) { + sim[k,,]<-sim[k,,]-prdMean+obsMean + } + } + if (method == "scaling") { + obsMean <- apply(obs, FUN = mean, MARGIN = c(2,3), na.rm = TRUE) + prdMean <- apply(pred, FUN = mean, MARGIN = c(2,3), na.rm = TRUE) + for (k in 1:dim(sim)[1]) { + sim[k,,]<-(sim[k,,]/prdMean)*obsMean + } + } + if (method == "qqmap") { + if (!any(grepl(varcode,c("pr","tp","precipitation","precip")))){ + for (i in 1:dim(sim)[2]) { + for (j in 1:dim(sim)[3]) { + if (any(!is.na(pred[,i,j])) & any(!is.na(obs[,i,j]))) { + ePrd <- ecdf(pred[,i,j]) + ################################################################################################3 + if(extrapolate == "constant"){ + exupsim <- which(sim[,i,j] > max(range(pred[,i,j], na.rm = TRUE))) + exdwnsim <- which(sim[,i,j] < min(range(pred[,i,j], na.rm = TRUE))) + ex <- c(exupsim,exdwnsim) + if (length(exupsim)>0){ + extmin <- max(pred[,i,j], na.rm = TRUE) + dif <- extmin - max(obs[,i,j], na.rm = TRUE) + sim[exupsim,i,j] <- sim[exupsim,i,j] - dif + } + if (length(exdwnsim)>0){ + extmin <- min(pred[,i,j], na.rm = TRUE) + dif <- extmin - min(obs[,i,j], na.rm = TRUE) + sim[exdwnsim,i,j] <- sim[exdwnsim,i,j] - dif + } + sim[-ex,i,j] <- quantile(obs[,i,j], probs = ePrd(sim[-ex,i,j]), na.rm = TRUE, type = 4) + }else{ + sim[,i,j] <- quantile(obs[,i,j], probs = ePrd(sim[,i,j]), na.rm = TRUE, type = 4) } + + ################################################################################################3 + } + } } - if (method == "qqmap") { - if (!any(grepl(varcode,c("pr","tp","precipitation","precip")))){ - for (i in 1:dim(sim)[2]) { - for (j in 1:dim(sim)[3]) { - if (any(!is.na(pred[,i,j])) & any(!is.na(obs[,i,j]))) { - ePrd<-ecdf(pred[,i,j]) - sim[,i,j]<-quantile(obs[,i,j], probs = ePrd(sim[,i,j]), na.rm = TRUE, type = 4) - } - } + } else { + for (i in 1:dim(sim)[2]) { + for (j in 1:dim(sim)[3]) { + if (any(!is.na(pred[,i,j])) & any(!is.na(obs[,i,j]))){ + if (length(which(pred[,i,j]>Pth[i,j]))>0){ + ## + ePrd<-ecdf(pred[which(pred[,i,j]>Pth[i,j]),i,j]) + ## + noRain<-which(sim[,i,j]<=Pth[i,j] & !is.na(sim[,i,j])) + rain<-which(sim[,i,j]>Pth[i,j] & !is.na(sim[,i,j])) + drizzle<-which(sim[,i,j]>Pth[i,j] & sim[,i,j] <= min(pred[which(pred[,i,j]>Pth[i,j]),i,j], na.rm = TRUE) & !is.na(sim[,i,j])) + if (length(rain)>0){ + eFrc<-ecdf(sim[rain,i,j]) #eFrc?? + + ##############################################################################3 + if(extrapolate == "constant"){ + exupsim <- which(sim[rain,i,j] > max(range(pred[,i,j], na.rm = TRUE))) # + exdwnsim <- which(sim[rain,i,j] < min(range(pred[,i,j], na.rm = TRUE))) + ex <- c(exupsim,exdwnsim) + if (length(exupsim)>0){ + extmin <- max(pred[,i,j], na.rm = TRUE) + dif <- extmin - max(obs[which(obs[,i,j] > threshold & !is.na(obs[,i,j])),i,j], na.rm = TRUE) + sim[rain[exupsim],i,j] <- sim[rain[exupsim],i,j] - dif } - } else { - for (i in 1:dim(sim)[2]) { - for (j in 1:dim(sim)[3]) { - if (any(!is.na(pred[,i,j])) & any(!is.na(obs[,i,j]))){ - if (length(which(pred[,i,j]>Pth[i,j]))>0){ - ePrd<-ecdf(pred[which(pred[,i,j]>Pth[i,j]),i,j]) - noRain<-which(sim[,i,j]<=Pth[i,j] & !is.na(sim[,i,j])) - rain<-which(sim[,i,j]>Pth[i,j] & !is.na(sim[,i,j])) - drizzle<-which(sim[,i,j]>Pth[i,j] & sim[,i,j]<=min(pred[which(pred[,i,j]>Pth[i,j]),i,j], na.rm = TRUE) & !is.na(sim[,i,j])) - if (length(rain)>0){ - eFrc<-ecdf(sim[rain,i,j]) - sim[rain,i,j]<-quantile(obs[which(obs[,i,j]>threshold & !is.na(obs[,i,j])),i,j], probs = ePrd(sim[rain,i,j]), na.rm = TRUE, type = 4) - } - if (length(drizzle)>0){ - sim[drizzle,i,j]<-quantile(sim[which(sim[,i,j]>min(pred[which(pred[,i,j]>Pth[i,j]),i,j], na.rm = TRUE) & !is.na(sim[,i,j])),i,j], probs = eFrc(sim[drizzle,i,j]), na.rm = TRUE, type = 4) - } - sim[noRain,i,j]<-0 - }else{ - noRain<-which(sim[,i,j]<=Pth[i,j] & !is.na(sim[,i,j])) - rain<-which(sim[,i,j]>Pth[i,j] & !is.na(sim[,i,j])) - if (length(rain)>0){ - eFrc<-ecdf(sim[rain,i,j]) - sim[rain,i,j]<-quantile(obs[which(obs[,i,j]>threshold & !is.na(obs[,i,j])),i,j], probs = ePrd(sim[rain,i,j]), na.rm = TRUE, type = 4) - } - sim[noRain,i,j]<-0 - } - } - } + if (length(exdwnsim)>0){ + extmin <- min(pred[,i,j], na.rm = TRUE) + dif <- extmin - min(obs[which(obs[,i,j] > threshold & !is.na(obs[,i,j])),i,j], na.rm = TRUE) + sim[rain[exdwnsim],i,j] <- sim[rain[exdwnsim],i,j] - dif } + sim[rain[-ex],i,j] <- quantile(obs[which(obs[,i,j] > threshold & !is.na(obs[,i,j])),i,j], probs = ePrd(sim[rain[-ex],i,j]), na.rm = TRUE, type = 4) + }else{ + sim[rain,i,j]<-quantile(obs[which(obs[,i,j]>threshold & !is.na(obs[,i,j])),i,j], probs = ePrd(sim[rain,i,j]), na.rm = TRUE, type = 4) + } + ################################################################################################3 + } + + if (length(drizzle)>0){ + + + sim[drizzle,i,j]<-quantile(sim[which(sim[,i,j]>min(pred[which(pred[,i,j]>Pth[i,j]),i,j], na.rm = TRUE) & !is.na(sim[,i,j])),i,j], probs = eFrc(sim[drizzle,i,j]), na.rm = TRUE, type = 4) + } + + sim[noRain,i,j]<-0 + }else{ + noRain<-which(sim[,i,j]<=Pth[i,j] & !is.na(sim[,i,j])) + rain<-which(sim[,i,j]>Pth[i,j] & !is.na(sim[,i,j])) + if (length(rain)>0){ + eFrc<-ecdf(sim[rain,i,j]) + sim[rain,i,j]<-quantile(obs[which(obs[,i,j] > threshold & !is.na(obs[,i,j])),i,j], probs = eFrc(sim[rain,i,j]), na.rm = TRUE, type = 4) + } + sim[noRain,i,j]<-0 } + } + } } - if (method == "piani" & any(grepl(varcode,c("pr","tp","precipitation","precip")))) { - for (i in 1:dim(sim)[2]){ - for (j in 1:dim(sim)[3]){ - if (nP[i,j]>0){ - ind<-which(obs[,i,j]>threshold & !is.na(obs[,i,j])) - obsGamma<-fitdistr(obs[ind,i,j],"gamma") - ind<-which(pred[,i,j]>0 & !is.na(pred[,i,j])) - prdGamma<-fitdistr(pred[ind,i,j],"gamma") - rain<-which(sim[,i,j]>Pth[i,j] & !is.na(sim[,i,j])) - noRain<-which(sim[,i,j]<=Pth[i,j] & !is.na(sim[,i,j])) - auxF<-pgamma(sim[rain,i,j],prdGamma$estimate[1], rate = prdGamma$estimate[2]) - sim[rain,i,j]<-qgamma(auxF,obsGamma$estimate[1], rate = obsGamma$estimate[2]) - sim[noRain,i,j]<-0 - } - } - } + } + } + if (method == "piani" & any(grepl(varcode,c("pr","tp","precipitation","precip")))) { + for (i in 1:dim(sim)[2]){ + for (j in 1:dim(sim)[3]){ + if (nP[i,j]>0){ + ind<-which(obs[,i,j]>threshold & !is.na(obs[,i,j])) + obsGamma<-fitdistr(obs[ind,i,j],"gamma") + ind<-which(pred[,i,j]>0 & !is.na(pred[,i,j])) + prdGamma<-fitdistr(pred[ind,i,j],"gamma") + rain<-which(sim[,i,j]>Pth[i,j] & !is.na(sim[,i,j])) + noRain<-which(sim[,i,j]<=Pth[i,j] & !is.na(sim[,i,j])) + auxF<-pgamma(sim[rain,i,j],prdGamma$estimate[1], rate = prdGamma$estimate[2]) + sim[rain,i,j]<-qgamma(auxF,obsGamma$estimate[1], rate = obsGamma$estimate[2]) + sim[noRain,i,j]<-0 + } + } + } + } + if (method == "gqm" & any(grepl(varcode,c("pr","tp","precipitation","precip")))) { + for (i in 1:dim(sim)[2]){ + for (j in 1:dim(sim)[3]){ + if (nP[i,j]>0){ + + ind<-which(obs[,i,j]>threshold & !is.na(obs[,i,j])) + + indgamma <- ind[which(obs[ind,i,j] < quantile(obs[ind,i,j], 0.95))] + indpareto <- ind[which(obs[ind,i,j] >= quantile(obs[ind,i,j], 0.95))] + + obsGQM <- fitdistr(obs[indgamma,i,j],"gamma") + obsGQM2 <- fpot(obs[indpareto,i,j], quantile(obs[ind,i,j], 0.95), "gpd", std.err = FALSE) + + ind <- which(pred[,i,j] > 0 & !is.na(pred[,i,j])) + indgamma <- ind[which(pred[ind,i,j]=quantile(pred[ind,i,j], 0.95))] + + prdGQM <- fitdistr(pred[indgamma,i,j], "gamma") + prdGQM2 <- fpot(pred[indpareto,i,j], quantile(pred[ind,i,j], 0.95), "gpd", std.err = FALSE) + + rain <- which(sim[,i,j] > Pth[i,j] & !is.na(sim[,i,j])) + noRain<-which(sim[,i,j] <= Pth[i,j] & !is.na(sim[,i,j])) + indgammasim <- rain[which(sim[rain,i,j] < quantile(pred[ind,i,j], 0.95))] + indparetosim <- rain[which(sim[rain,i,j] >= quantile(pred[ind,i,j], 0.95))] + + auxF <- pgamma(sim[indgammasim,i,j],prdGQM$estimate[1], rate = prdGQM$estimate[2]) + auxF2 <- pgpd(sim[indparetosim,i,j],loc = 0, scale = prdGQM2$estimate[1], shape = prdGQM2$estimate[2]) + + sim[indgammasim,i,j] <- qgamma(auxF, obsGQM$estimate[1], rate = obsGQM$estimate[2]) + sim[indparetosim[which(auxF2<1)],i,j] <- qgpd(auxF2[which(auxF2 < 1)], loc = 0, scale = obsGQM2$estimate[1], shape = obsGQM2$estimate[2]) + sim[indparetosim[which(auxF2==1)],i,j] <- max(obs[indpareto,i,j], na.rm = TRUE) + + sim[noRain,i,j]<-0 + } } - return(sim) + } + } + return(sim) } + # End From d806a374fdb4c715fd490fa0192fa7c10e57f08a Mon Sep 17 00:00:00 2001 From: jbedia Date: Fri, 26 Jun 2015 13:47:40 +0200 Subject: [PATCH 3/5] Bug fix in subdaily data def Whole-year seasons (season = 1:12) are now included in a dateSliceList structure like ordinary seasons for compatibility --- R/getTimeDomain.R | 33 +++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/R/getTimeDomain.R b/R/getTimeDomain.R index 1ada9d76..7d04ffcf 100644 --- a/R/getTimeDomain.R +++ b/R/getTimeDomain.R @@ -86,27 +86,32 @@ getTimeDomain <- function(grid, dic, season, years, time, aggr.d, aggr.m) { } else { timeInd <- which((timeDates$year + 1900) %in% years & (timeDates$mon + 1) %in% season) } - dateSlice <- timeDates[timeInd] + timeDates <- timeDates[timeInd] timeIndList <- list() dateSliceList <- list() - if (length(dateSlice) > 1) { + if (length(timeDates) > 1) { brkInd <- rep(1, length(timeInd)) for (i in 2:length(timeInd)) { brkInd[i] <- timeInd[i] - timeInd[i-1] } - brkInd <- c(1, which(brkInd > 1), length(timeInd) + 1) - if (length(brkInd) == 0) { - timeIndList[[1]] <- timeInd - 1 - dateSliceList[[1]] <- dateSlice - } else { - for (i in 2:length(brkInd)) { - timeIndList[[i - 1]] <- timeInd[brkInd[i - 1] : (brkInd[i] - 1)] - 1 - dateSliceList[[i - 1]] <- dateSlice[brkInd[i - 1] : (brkInd[i] - 1)] - } - } - } else { + if (length(which(brkInd > 1)) != 0) { + brkInd <- c(1, which(brkInd > 1), length(timeInd) + 1) + if (length(brkInd) == 0) { + timeIndList[[1]] <- timeInd - 1 + dateSliceList[[1]] <- timeDates + } else { + for (i in 2:length(brkInd)) { + timeIndList[[i - 1]] <- timeInd[brkInd[i - 1] : (brkInd[i] - 1)] - 1 + dateSliceList[[i - 1]] <- timeDates[brkInd[i - 1] : (brkInd[i] - 1)] + } + } + } else { + dateSliceList <- lapply(unique(timeDates$year), function(x) timeDates[which(timeDates$year == x)]) + timeIndList <- lapply(unique(timeDates$year), function(x) timeInd[which(timeDates$year == x)]) + } + } else { timeIndList[[1]] <- timeInd - 1 - dateSliceList[[1]] <- dateSlice + dateSliceList[[1]] <- timeDates } if (time == "DD" | time == "none") { timeStride <- 1L From 3b8ef6795f7dcdc1519a9b4a573ab569223a374f Mon Sep 17 00:00:00 2001 From: jbedia Date: Fri, 26 Jun 2015 14:14:08 +0200 Subject: [PATCH 4/5] Update previous to minor v release --- DESCRIPTION | 14 +++++++------ NAMESPACE | 3 +++ R/biasCorrection.R | 32 +++++++++++++++++++++------- man/biasCorrection.Rd | 49 ++++++++++++++++++++++++++++++------------- man/calibrateProj.Rd | 13 ++++++------ 5 files changed, 78 insertions(+), 33 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index fc85a53a..2e691ab1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -12,26 +12,28 @@ Imports: MASS, verification, scales, + evd Suggests: ecomsUDG.Raccess, Type: Package Title: Climate data manipulation and statistical downscaling -Version: 0.7-0 -Date: 15-May-2015 +Version: 0.8-0 +Date: 26-Jun-2015 Authors@R: as.person(c( "Joaquin Bedia [ctb, cre]", + "Antonio Cofino [ctb]", "Sixto Herrera [ctb]", - "Maria Dolores Frias [ctb]", "Jesus Fernandez [ctb]", "Wietse Franssen [ctb]", + "Maria Dolores Frias [ctb]", + "Maialen Iturbide [ctb]", "Max Tuni [ctb]", - "Antonio Cofino [ctb]", "Santander Meteorology Group [aut]")) BugReports: https://github.com/SantanderMetGroup/downscaleR/issues URL: https://github.com/SantanderMetGroup/downscaleR/wiki Description: Load climate and weather data into R and performs climate data analysis and visualization, including model calibration (bias correction, - qq-mapping...), and perfect-prog statistical downscaling. The package is - conceived for dealing also with forecast (multi-member) data. + qq-mapping...), and perfect-prog statistical downscaling. Focused on daily data, + the package is conceived for dealing also with forecast (multi-member) data. License: GPL (>= 3) LazyData: true diff --git a/NAMESPACE b/NAMESPACE index a4958d69..f1ddfdf5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -50,6 +50,9 @@ importFrom(abind,abind) importFrom(abind,asub) importFrom(downscaleR.java,javaCalendarDate2rPOSIXlt) importFrom(downscaleR.java,javaString2rChar) +importFrom(evd,fpot) +importFrom(evd,pgpd) +importFrom(evd,qgpd) importFrom(fields,image.plot) importFrom(fields,interp.surface.grid) importFrom(fields,rdist) diff --git a/R/biasCorrection.R b/R/biasCorrection.R index 5bfd344b..7a2c3cfd 100644 --- a/R/biasCorrection.R +++ b/R/biasCorrection.R @@ -80,17 +80,33 @@ #' @export #' @examples \dontrun{ #' # These are the paths to the package built-in GSN and NCEP datasets -#' gsn.data.dir <- file.path(find.package("downscaleR"), "datasets/observations/GSN_Iberia") -#' ncep.data.dir <- file.path(find.package("downscaleR"), "datasets/reanalysis/Iberia_NCEP/Iberia_NCEP.ncml") +#' gsn.data.dir <- file.path(find.package("downscaleR.java"), +#' "datasets/observations/GSN_Iberia") +#' ncep.data.dir <- file.path(find.package("downscaleR.java"), +#' "datasets/reanalysis/Iberia_NCEP/Iberia_NCEP.ncml") #' gsn.inv <- dataInventory(gsn.data.dir) #' ncep.inv <- dataInventory(ncep.data.dir) #' str(gsn.inv) #' str(ncep.inv) -#' # Load precipitation for boreal winter (DJF) in the train (1991-2000) and test (2001-2010) periods, for the observations (GSN_Iberia) and the Iberia_NCEP datasets -#' obs <- loadStationData(dataset = gsn.data.dir, var="precip", lonLim = c(-12,10), latLim = c(33,47), season=c(12,1,2), years = 1991:2000) -#' prd <- loadGridData(ncep.data.dir, var = "tp", lonLim = c(-12,10), latLim = c(33,47), season = c(12,1,2), years = 1991:2000) -#' sim <- loadGridData(ncep.data.dir, var = "tp", lonLim = c(-12,10), latLim = c(33,47), season = c(12,1,2), years = 2001:2010) -#' # Interpolation of the observations onto the grid of model: we use the method "nearest" and the getGrid function to ensure spatial consistency: +#' # Load precipitation for boreal winter (DJF) in the train (1991-2000) and test (2001-2010) periods, +#' # for the observations (GSN_Iberia) and the Iberia_NCEP datasets +#' obs <- loadStationData(dataset = gsn.data.dir, +#' var="precip", +#' lonLim = c(-12,10), latLim = c(33,47), +#' season = c(12,1,2), +#' years = 1991:2000) +#' prd <- loadGridData(dataset = ncep.data.dir, +#' var = "tp", +#' lonLim = c(-12,10), latLim = c(33,47), +#' season = c(12,1,2), +#' years = 1991:2000) +#' sim <- loadGridData(dataset = ncep.data.dir, +#' var = "tp", +#' lonLim = c(-12,10), latLim = c(33,47), +#' season = c(12,1,2), +#' years = 2001:2010) +#' # Interpolation of the observations onto the grid of model: we use the method "nearest" and the +#' # 'getGrid' function to ensure spatial consistency: #' obs <- interpGridData(obs, new.grid = getGrid(prd), method = "nearest") #' # Apply the bias correction method: #' simBC <- biasCorrection (obs, prd, sim, method = "qqmap", pr.threshold = 1) # qq-mapping @@ -304,6 +320,8 @@ biasCorrection <- function (obs, pred, sim, method = c("qqmap", "delta", "scalin #' @family calibration #' @importFrom MASS fitdistr #' @importFrom evd fpot +#' @importFrom evd qgpd +#' @importFrom evd pgpd #' @keywords internal #' diff --git a/man/biasCorrection.Rd b/man/biasCorrection.Rd index c8891c9f..9a91c5e1 100644 --- a/man/biasCorrection.Rd +++ b/man/biasCorrection.Rd @@ -5,8 +5,8 @@ \title{Bias correction methods} \usage{ biasCorrection(obs, pred, sim, method = c("qqmap", "delta", "scaling", - "unbiasing", "piani"), pr.threshold = 1, multi.member = TRUE, - window = NULL) + "unbiasing", "piani", "gqm"), pr.threshold = 1, multi.member = TRUE, + window = NULL, extrapolation = c("no", "constant")) } \arguments{ \item{obs}{A field or station data containing the observed climate data for the training period} @@ -19,7 +19,7 @@ the same variable as \code{obs}, in the case of model calibration (bias correcti \item{sim}{A field containing the simulated climate for the variables used in \code{pred}, but considering the test period.} \item{method}{method applied. Current accepted values are \code{"qqmap"}, \code{"delta"}, -\code{"scaling"}, \code{"unbiasing"} and \code{"piani"}. See details.} +\code{"scaling"}, \code{"unbiasing"}, \code{"piani"} and \code{"gqm"}. See details.} \item{pr.threshold}{The minimum value that is considered as a non-zero precipitation. Ignored for \code{varcode} values different from \code{"pr"}. Default to 1 (assuming mm).} @@ -29,6 +29,10 @@ Ignored if the dataset has no members.} \item{window}{Numeric value specifying the time window width used to calibrate. The window is centered on the target day. Default to \code{NULL}, which considers the whole period available.} + +\item{extrapolation}{Character indicating the extrapolation method to be applied to correct values in +\code{"sim"} that are out of the range of \code{"pred"}. Extrapolation is applied only to the \code{"qqmap"} method, +thus, this argument is ignored if other bias correction method is selected. Default is \code{"no"} (do not extrapolate).} } \value{ A calibrated object of the same spatio-temporal extent of the input field @@ -38,7 +42,7 @@ Implementation of several standard bias correction methods } \details{ The methods available are \code{"qqmap"}, \code{"delta"}, \code{"unbiasing"}, -\code{"scaling"} and \code{"Piani"} (the latter used only for precipitation). +\code{"scaling"}, \code{"Piani"}, \code{"gqm"} (the two latter used only for precipitation). Next we make a brief description of each method: \strong{Delta} @@ -69,28 +73,43 @@ This method is applicable to any kind of variable. This method is described in Piani et al. 2010 and is applicable only to precipitation. It is based on the initial assumption that both observed and simulated intensity distributions are well approximated by the gamma distribution, therefore is a parametric q-q map that uses the theorical instead of the empirical distribution. + +\strong{Generalized Quantile Mapping (gqm)} + +This method is described in Gutjahr and Heinemann 2013. It is applicable only to precipitation and is similar to the Piani method. It applies a +gamma distribution to values under the threshold given by the 95th percentile (proosed by Yang et al. 2010) and a general Pareto +distribution (GPD) to values above the threshold. } \examples{ \dontrun{ # These are the paths to the package built-in GSN and NCEP datasets gsn.data.dir <- file.path(find.package("downscaleR.java"), - "datasets/observations/GSN_Iberia") + "datasets/observations/GSN_Iberia") ncep.data.dir <- file.path(find.package("downscaleR.java"), - "datasets/reanalysis/Iberia_NCEP/Iberia_NCEP.ncml") + "datasets/reanalysis/Iberia_NCEP/Iberia_NCEP.ncml") gsn.inv <- dataInventory(gsn.data.dir) ncep.inv <- dataInventory(ncep.data.dir) str(gsn.inv) str(ncep.inv) # Load precipitation for boreal winter (DJF) in the train (1991-2000) and test (2001-2010) periods, # for the observations (GSN_Iberia) and the Iberia_NCEP datasets -obs <- loadStationData(dataset = gsn.data.dir, var="precip", lonLim = c(-12,10), latLim = c(33,47), - season=c(12,1,2), years = 1991:2000) -prd <- loadGridData(ncep.data.dir, var = "tp", lonLim = c(-12,10), latLim = c(33,47), - season = c(12,1,2), years = 1991:2000) -sim <- loadGridData(ncep.data.dir, var = "tp", lonLim = c(-12,10), latLim = c(33,47), - season = c(12,1,2), years = 2001:2010) -# Interpolation of the observations onto the grid of model: we use the method "nearest" -# and the getGrid function to ensure spatial consistency: +obs <- loadStationData(dataset = gsn.data.dir, + var="precip", + lonLim = c(-12,10), latLim = c(33,47), + season = c(12,1,2), + years = 1991:2000) +prd <- loadGridData(dataset = ncep.data.dir, + var = "tp", + lonLim = c(-12,10), latLim = c(33,47), + season = c(12,1,2), + years = 1991:2000) +sim <- loadGridData(dataset = ncep.data.dir, + var = "tp", + lonLim = c(-12,10), latLim = c(33,47), + season = c(12,1,2), + years = 2001:2010) +# Interpolation of the observations onto the grid of model: we use the method "nearest" and the +# 'getGrid' function to ensure spatial consistency: obs <- interpGridData(obs, new.grid = getGrid(prd), method = "nearest") # Apply the bias correction method: simBC <- biasCorrection (obs, prd, sim, method = "qqmap", pr.threshold = 1) # qq-mapping @@ -109,6 +128,8 @@ S. Herrera \email{sixto@predictia.es} \item A. Amengual, V. Homar, R. Romero, S. Alonso, and C. Ramis (2012) A Statistical Adjustment of Regional Climate Model Outputs to Local Scales: Application to Platja de Palma, Spain. J. Clim., 25, 939-957 \item C. Piani, J. O. Haerter and E. Coppola (2009) Statistical bias correction for daily precipitation in regional climate models over Europe, Theoretical and Applied Climatology, 99, 187-192 + +\item O. Gutjahr and G. Heinemann (2013) Comparing precipitation bias correction methods for high-resolution regional climate simulations using COSMO-CLM, Theoretical and Applied Climatology, 114, 511-529 } } \seealso{ diff --git a/man/calibrateProj.Rd b/man/calibrateProj.Rd index 1f925213..e3979b75 100644 --- a/man/calibrateProj.Rd +++ b/man/calibrateProj.Rd @@ -5,8 +5,8 @@ \title{Bias correction methods} \usage{ calibrateProj(obs, pred, sim, method = c("qqmap", "delta", "scaling", - "unbiasing", "piani"), varcode = c("tas", "hurs", "tp", "pr", "wss"), - pr.threshold = 1) + "unbiasing", "piani", "gqm"), varcode = c("tas", "hurs", "tp", "pr", "wss"), + pr.threshold = 1, extrapolate = c("no", "constant")) } \arguments{ \item{obs}{A field or station data containing the observed climate data for the training period} @@ -19,7 +19,7 @@ the same variable as \code{obs}, in the case of model calibration (bias correcti \item{sim}{A field containing the simulated climate for the variables used in \code{pred}, but considering the test period.} \item{method}{method applied. Current accepted values are \code{"qqmap"}, \code{"delta"}, -\code{"scaling"}, \code{"unbiasing"} and \code{"piani"}. See details.} +\code{"scaling"}, \code{"unbiasing"}, \code{"piani"} and \code{"gqm"}. See details.} \item{varcode}{Variable code. This is not the variable itself to be corrected, but rather it referes to its nature and distributional properties. For instance, \code{"tas"} applies for @@ -30,13 +30,14 @@ but without an upper bound) and \code{"pr"}, specifically for precipitation.} \item{pr.threshold}{The minimum value that is considered as a non-zero precipitation. Ignored for \code{varcode} values different from \code{"pr"}. Default to 1 (assuming mm).} + +\item{extrapolate}{Character indicating the extrapolation method to be applied to correct values in +\code{"sim"} that are out of the range of \code{"pred"}. Extrapolation is applied only to the \code{"qqmap"} method, +thus, this argument is ignored if other bias correction method is selected. Default is \code{"no"} (do not extrapolate).} } \description{ Implementation of several standard bias correction methods } -\author{ -S. Herrera \email{sixto@predictia.es} -} \seealso{ Other downscaling: \code{\link{analogs}}; \code{\link{biasCorrection}}; \code{\link{glimpr}}; From 97f8eb8f429be2091b2581b0e0bbef654fc63cc8 Mon Sep 17 00:00:00 2001 From: jbedia Date: Fri, 26 Jun 2015 14:14:08 +0200 Subject: [PATCH 5/5] Update previous to minor v release --- DESCRIPTION | 14 +++++++------ NAMESPACE | 3 +++ NEWS | 10 +++------ R/biasCorrection.R | 32 +++++++++++++++++++++------- man/biasCorrection.Rd | 49 ++++++++++++++++++++++++++++++------------- man/calibrateProj.Rd | 13 ++++++------ 6 files changed, 81 insertions(+), 40 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index fc85a53a..2e691ab1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -12,26 +12,28 @@ Imports: MASS, verification, scales, + evd Suggests: ecomsUDG.Raccess, Type: Package Title: Climate data manipulation and statistical downscaling -Version: 0.7-0 -Date: 15-May-2015 +Version: 0.8-0 +Date: 26-Jun-2015 Authors@R: as.person(c( "Joaquin Bedia [ctb, cre]", + "Antonio Cofino [ctb]", "Sixto Herrera [ctb]", - "Maria Dolores Frias [ctb]", "Jesus Fernandez [ctb]", "Wietse Franssen [ctb]", + "Maria Dolores Frias [ctb]", + "Maialen Iturbide [ctb]", "Max Tuni [ctb]", - "Antonio Cofino [ctb]", "Santander Meteorology Group [aut]")) BugReports: https://github.com/SantanderMetGroup/downscaleR/issues URL: https://github.com/SantanderMetGroup/downscaleR/wiki Description: Load climate and weather data into R and performs climate data analysis and visualization, including model calibration (bias correction, - qq-mapping...), and perfect-prog statistical downscaling. The package is - conceived for dealing also with forecast (multi-member) data. + qq-mapping...), and perfect-prog statistical downscaling. Focused on daily data, + the package is conceived for dealing also with forecast (multi-member) data. License: GPL (>= 3) LazyData: true diff --git a/NAMESPACE b/NAMESPACE index a4958d69..f1ddfdf5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -50,6 +50,9 @@ importFrom(abind,abind) importFrom(abind,asub) importFrom(downscaleR.java,javaCalendarDate2rPOSIXlt) importFrom(downscaleR.java,javaString2rChar) +importFrom(evd,fpot) +importFrom(evd,pgpd) +importFrom(evd,qgpd) importFrom(fields,image.plot) importFrom(fields,interp.surface.grid) importFrom(fields,rdist) diff --git a/NEWS b/NEWS index 69cff356..8dfadcb5 100644 --- a/NEWS +++ b/NEWS @@ -1,11 +1,7 @@ -downscaleR 0.7-0 +downscaleR 0.8-0 ================ -* Improved time aggregation capabilities - * Added on-the-fly monthly aggregation (time = "MM") - * Aggregation functions can be indicated for both daily and mnonthly aggregations -* netCDF-4 export removed from package -* Built-in datasets removed and moved to "downscaleR.java" package -* New variable metadata included (units, temporal aggregation info...) +* New Generalized Quantile Mapping method (Gutjahr and Heinemann 2013) +* New extrapolation feature of Quantile-quantile mapping method for unprecedent values in the simulation period * Minor bug fixes and documentation improvements diff --git a/R/biasCorrection.R b/R/biasCorrection.R index 5bfd344b..7a2c3cfd 100644 --- a/R/biasCorrection.R +++ b/R/biasCorrection.R @@ -80,17 +80,33 @@ #' @export #' @examples \dontrun{ #' # These are the paths to the package built-in GSN and NCEP datasets -#' gsn.data.dir <- file.path(find.package("downscaleR"), "datasets/observations/GSN_Iberia") -#' ncep.data.dir <- file.path(find.package("downscaleR"), "datasets/reanalysis/Iberia_NCEP/Iberia_NCEP.ncml") +#' gsn.data.dir <- file.path(find.package("downscaleR.java"), +#' "datasets/observations/GSN_Iberia") +#' ncep.data.dir <- file.path(find.package("downscaleR.java"), +#' "datasets/reanalysis/Iberia_NCEP/Iberia_NCEP.ncml") #' gsn.inv <- dataInventory(gsn.data.dir) #' ncep.inv <- dataInventory(ncep.data.dir) #' str(gsn.inv) #' str(ncep.inv) -#' # Load precipitation for boreal winter (DJF) in the train (1991-2000) and test (2001-2010) periods, for the observations (GSN_Iberia) and the Iberia_NCEP datasets -#' obs <- loadStationData(dataset = gsn.data.dir, var="precip", lonLim = c(-12,10), latLim = c(33,47), season=c(12,1,2), years = 1991:2000) -#' prd <- loadGridData(ncep.data.dir, var = "tp", lonLim = c(-12,10), latLim = c(33,47), season = c(12,1,2), years = 1991:2000) -#' sim <- loadGridData(ncep.data.dir, var = "tp", lonLim = c(-12,10), latLim = c(33,47), season = c(12,1,2), years = 2001:2010) -#' # Interpolation of the observations onto the grid of model: we use the method "nearest" and the getGrid function to ensure spatial consistency: +#' # Load precipitation for boreal winter (DJF) in the train (1991-2000) and test (2001-2010) periods, +#' # for the observations (GSN_Iberia) and the Iberia_NCEP datasets +#' obs <- loadStationData(dataset = gsn.data.dir, +#' var="precip", +#' lonLim = c(-12,10), latLim = c(33,47), +#' season = c(12,1,2), +#' years = 1991:2000) +#' prd <- loadGridData(dataset = ncep.data.dir, +#' var = "tp", +#' lonLim = c(-12,10), latLim = c(33,47), +#' season = c(12,1,2), +#' years = 1991:2000) +#' sim <- loadGridData(dataset = ncep.data.dir, +#' var = "tp", +#' lonLim = c(-12,10), latLim = c(33,47), +#' season = c(12,1,2), +#' years = 2001:2010) +#' # Interpolation of the observations onto the grid of model: we use the method "nearest" and the +#' # 'getGrid' function to ensure spatial consistency: #' obs <- interpGridData(obs, new.grid = getGrid(prd), method = "nearest") #' # Apply the bias correction method: #' simBC <- biasCorrection (obs, prd, sim, method = "qqmap", pr.threshold = 1) # qq-mapping @@ -304,6 +320,8 @@ biasCorrection <- function (obs, pred, sim, method = c("qqmap", "delta", "scalin #' @family calibration #' @importFrom MASS fitdistr #' @importFrom evd fpot +#' @importFrom evd qgpd +#' @importFrom evd pgpd #' @keywords internal #' diff --git a/man/biasCorrection.Rd b/man/biasCorrection.Rd index c8891c9f..9a91c5e1 100644 --- a/man/biasCorrection.Rd +++ b/man/biasCorrection.Rd @@ -5,8 +5,8 @@ \title{Bias correction methods} \usage{ biasCorrection(obs, pred, sim, method = c("qqmap", "delta", "scaling", - "unbiasing", "piani"), pr.threshold = 1, multi.member = TRUE, - window = NULL) + "unbiasing", "piani", "gqm"), pr.threshold = 1, multi.member = TRUE, + window = NULL, extrapolation = c("no", "constant")) } \arguments{ \item{obs}{A field or station data containing the observed climate data for the training period} @@ -19,7 +19,7 @@ the same variable as \code{obs}, in the case of model calibration (bias correcti \item{sim}{A field containing the simulated climate for the variables used in \code{pred}, but considering the test period.} \item{method}{method applied. Current accepted values are \code{"qqmap"}, \code{"delta"}, -\code{"scaling"}, \code{"unbiasing"} and \code{"piani"}. See details.} +\code{"scaling"}, \code{"unbiasing"}, \code{"piani"} and \code{"gqm"}. See details.} \item{pr.threshold}{The minimum value that is considered as a non-zero precipitation. Ignored for \code{varcode} values different from \code{"pr"}. Default to 1 (assuming mm).} @@ -29,6 +29,10 @@ Ignored if the dataset has no members.} \item{window}{Numeric value specifying the time window width used to calibrate. The window is centered on the target day. Default to \code{NULL}, which considers the whole period available.} + +\item{extrapolation}{Character indicating the extrapolation method to be applied to correct values in +\code{"sim"} that are out of the range of \code{"pred"}. Extrapolation is applied only to the \code{"qqmap"} method, +thus, this argument is ignored if other bias correction method is selected. Default is \code{"no"} (do not extrapolate).} } \value{ A calibrated object of the same spatio-temporal extent of the input field @@ -38,7 +42,7 @@ Implementation of several standard bias correction methods } \details{ The methods available are \code{"qqmap"}, \code{"delta"}, \code{"unbiasing"}, -\code{"scaling"} and \code{"Piani"} (the latter used only for precipitation). +\code{"scaling"}, \code{"Piani"}, \code{"gqm"} (the two latter used only for precipitation). Next we make a brief description of each method: \strong{Delta} @@ -69,28 +73,43 @@ This method is applicable to any kind of variable. This method is described in Piani et al. 2010 and is applicable only to precipitation. It is based on the initial assumption that both observed and simulated intensity distributions are well approximated by the gamma distribution, therefore is a parametric q-q map that uses the theorical instead of the empirical distribution. + +\strong{Generalized Quantile Mapping (gqm)} + +This method is described in Gutjahr and Heinemann 2013. It is applicable only to precipitation and is similar to the Piani method. It applies a +gamma distribution to values under the threshold given by the 95th percentile (proosed by Yang et al. 2010) and a general Pareto +distribution (GPD) to values above the threshold. } \examples{ \dontrun{ # These are the paths to the package built-in GSN and NCEP datasets gsn.data.dir <- file.path(find.package("downscaleR.java"), - "datasets/observations/GSN_Iberia") + "datasets/observations/GSN_Iberia") ncep.data.dir <- file.path(find.package("downscaleR.java"), - "datasets/reanalysis/Iberia_NCEP/Iberia_NCEP.ncml") + "datasets/reanalysis/Iberia_NCEP/Iberia_NCEP.ncml") gsn.inv <- dataInventory(gsn.data.dir) ncep.inv <- dataInventory(ncep.data.dir) str(gsn.inv) str(ncep.inv) # Load precipitation for boreal winter (DJF) in the train (1991-2000) and test (2001-2010) periods, # for the observations (GSN_Iberia) and the Iberia_NCEP datasets -obs <- loadStationData(dataset = gsn.data.dir, var="precip", lonLim = c(-12,10), latLim = c(33,47), - season=c(12,1,2), years = 1991:2000) -prd <- loadGridData(ncep.data.dir, var = "tp", lonLim = c(-12,10), latLim = c(33,47), - season = c(12,1,2), years = 1991:2000) -sim <- loadGridData(ncep.data.dir, var = "tp", lonLim = c(-12,10), latLim = c(33,47), - season = c(12,1,2), years = 2001:2010) -# Interpolation of the observations onto the grid of model: we use the method "nearest" -# and the getGrid function to ensure spatial consistency: +obs <- loadStationData(dataset = gsn.data.dir, + var="precip", + lonLim = c(-12,10), latLim = c(33,47), + season = c(12,1,2), + years = 1991:2000) +prd <- loadGridData(dataset = ncep.data.dir, + var = "tp", + lonLim = c(-12,10), latLim = c(33,47), + season = c(12,1,2), + years = 1991:2000) +sim <- loadGridData(dataset = ncep.data.dir, + var = "tp", + lonLim = c(-12,10), latLim = c(33,47), + season = c(12,1,2), + years = 2001:2010) +# Interpolation of the observations onto the grid of model: we use the method "nearest" and the +# 'getGrid' function to ensure spatial consistency: obs <- interpGridData(obs, new.grid = getGrid(prd), method = "nearest") # Apply the bias correction method: simBC <- biasCorrection (obs, prd, sim, method = "qqmap", pr.threshold = 1) # qq-mapping @@ -109,6 +128,8 @@ S. Herrera \email{sixto@predictia.es} \item A. Amengual, V. Homar, R. Romero, S. Alonso, and C. Ramis (2012) A Statistical Adjustment of Regional Climate Model Outputs to Local Scales: Application to Platja de Palma, Spain. J. Clim., 25, 939-957 \item C. Piani, J. O. Haerter and E. Coppola (2009) Statistical bias correction for daily precipitation in regional climate models over Europe, Theoretical and Applied Climatology, 99, 187-192 + +\item O. Gutjahr and G. Heinemann (2013) Comparing precipitation bias correction methods for high-resolution regional climate simulations using COSMO-CLM, Theoretical and Applied Climatology, 114, 511-529 } } \seealso{ diff --git a/man/calibrateProj.Rd b/man/calibrateProj.Rd index 1f925213..e3979b75 100644 --- a/man/calibrateProj.Rd +++ b/man/calibrateProj.Rd @@ -5,8 +5,8 @@ \title{Bias correction methods} \usage{ calibrateProj(obs, pred, sim, method = c("qqmap", "delta", "scaling", - "unbiasing", "piani"), varcode = c("tas", "hurs", "tp", "pr", "wss"), - pr.threshold = 1) + "unbiasing", "piani", "gqm"), varcode = c("tas", "hurs", "tp", "pr", "wss"), + pr.threshold = 1, extrapolate = c("no", "constant")) } \arguments{ \item{obs}{A field or station data containing the observed climate data for the training period} @@ -19,7 +19,7 @@ the same variable as \code{obs}, in the case of model calibration (bias correcti \item{sim}{A field containing the simulated climate for the variables used in \code{pred}, but considering the test period.} \item{method}{method applied. Current accepted values are \code{"qqmap"}, \code{"delta"}, -\code{"scaling"}, \code{"unbiasing"} and \code{"piani"}. See details.} +\code{"scaling"}, \code{"unbiasing"}, \code{"piani"} and \code{"gqm"}. See details.} \item{varcode}{Variable code. This is not the variable itself to be corrected, but rather it referes to its nature and distributional properties. For instance, \code{"tas"} applies for @@ -30,13 +30,14 @@ but without an upper bound) and \code{"pr"}, specifically for precipitation.} \item{pr.threshold}{The minimum value that is considered as a non-zero precipitation. Ignored for \code{varcode} values different from \code{"pr"}. Default to 1 (assuming mm).} + +\item{extrapolate}{Character indicating the extrapolation method to be applied to correct values in +\code{"sim"} that are out of the range of \code{"pred"}. Extrapolation is applied only to the \code{"qqmap"} method, +thus, this argument is ignored if other bias correction method is selected. Default is \code{"no"} (do not extrapolate).} } \description{ Implementation of several standard bias correction methods } -\author{ -S. Herrera \email{sixto@predictia.es} -} \seealso{ Other downscaling: \code{\link{analogs}}; \code{\link{biasCorrection}}; \code{\link{glimpr}};