#' Biascorrect the input timeseries or hyfo dataset
#'
#' Biascorrect the input time series or dataset, the input time series or dataset should consist of observation, hindcast, and forecast.
#' observation and hindcast should belong to the same period, in order to calibrate. Then the modified forecast
#' will be returned. If the input is a time series, first column should be date column and rest columns should be
#' the value column. If the input is a hyfo dataset, the dataset should be the result of \code{loadNcdf}, or a list
#' file with the same format.
#'
#'
#' @param frc a hyfo grid data output or a dataframe (time series) consists of Date column and one or more value columns,
#' representing the forecast to be calibrated.
#' @param hindcast a hyfo grid data output or a dataframe(time series) consists of Date column and one or more value columns,
#' representing the hindcast data. This data will be used in the calibration of the forecast, so it's better to have the same date period as
#' observation data. Check details for more information.
#' @param obs a hyfo grid data output or a dataframe (time series) consists of Date column and one or more value columns,
#' representing the observation data.
#' @param method bias correct method, including 'delta', 'scaling'..., default is 'scaling'
#' @param scaleType only when the method "scaling" is chosen, scaleType will be available. Two different types
#' of scaling method, 'add' and 'multi', which means additive and multiplicative scaling method. More info check
#' details. Default scaleType is 'multi'.
#' @param preci If the precipitation is biascorrected, then you have to assign \code{preci = TRUE}. Since for
#' precipitation, some biascorrect methods may not apply to, or some methods are specially for precipitation.
#' Default is FALSE, refer to details.
#' @param prThreshold The minimum value that is considered as a non-zero precipitation. Default to 0 (assuming mm). If you want
#' to use precipitation biascorrect, you should consider carefully how to set this threshold, usually is 1. But you
#' can try with different numbers to see the results.
#' @param extrapolate When use 'eqm' method, and 'no' is set, modified frc is bounded by the range of obs.
#' If 'constant' is set, modified frc is not bounded by the range of obs. Default is 'no'.
#' @details
#'
#' Since climate forecast is based on global condition, when downscaling to different regions, it may include
#' some bias, biascorrection is used then to fix the bias.
#'
#' \strong{Hindcast}
#'
#' In order to bias correct, we need to pick up some data from the forecast to train with
#' the observation, which is called hindcast in this function. Using hindcast and observation,
#' the program can analyze the bias and correct the bias in the forecast.
#'
#' Hindcast should have \strong{EVERY} attributes that forecast has.
#'
#' Hindcast is also called re-forecast, is the forecast of the past. E.g. you have a forecast from year 2000-2010, assuming now you are in 2005. So from 2000-2005, this period
#' is the hindcast period, and 2005-2010, this period is the forecast period.
#'
#' Hindcast can be the same as forecast, i.e., you can use forecast itself as hindcast to train the bias correction.
#'
#'
#' \strong{How it works}
#'
#' Forecast product has to be calibrated, usually the system is doing forecast in real time. So, e.g., if the
#' forecast starts from year 2000, assuming you are in year 2003, then you will have 3 years' hindcast
#' data (year 2000-2003), which can be used to calibrate. And your forecast period is (2003-2004)
#'
#' E.g. you have observation from 2001-2002, this is your input obs. Then you can take the same
#' period (2001-2002) from the forecast, which is the hindcast period. For forecast, you can take any period.
#' The program will evaluate the obs and hindcast, to get the modification of the forecast, and then add the
#' modification to the forecast data.
#'
#' The more categorized input, the more accurate result you will get. E.g., if you want to
#' bias correct a forecast for winter season. So you'd better to extract all the winter period
#' in the hindcast and observation to train. \code{extractPeriod} can be used for this purpose.
#'
#' \strong{method}
#'
#' Different methods used in the bias correction. Among which, delta, scaling can be applied
#' to different kinds of parameters, with no need to set \code{preci}; eqm has two conditions for rainfall data and other data,
#' it needs user to input \code{preci = TRUE/FALSE} to point to different conditions; gqm is
#' designed for rainfall data, so \code{preci = TRUE} needs to be set.
#'
#' \strong{delta}
#'
#' This method consists on adding to the observations the mean change signal (delta method).
#' This method is applicable to any kind of variable but it is preferable to avoid it for bounded variables
#' (e.g. precipitation, wind speed, etc.) because values out of the variable range could be obtained
#' (e.g. negative wind speeds...)
#'
#' \strong{scaling}
#'
#' This method consists on scaling the simulation with the difference (additive) or quotient (multiplicative)
#' between the observed and simulated means in the train period. The \code{additive} or \code{multiplicative}
#' correction is defined by parameter \code{scaling.type} (default is \code{additive}).
#' The additive version is preferably applicable to unbounded variables (e.g. temperature)
#' and the multiplicative to variables with a lower bound (e.g. precipitation, because it also preserves the frequency).
#'
#' \strong{eqm}
#'
#' Empirical Quantile Mapping. This is a very extended bias correction method which consists on calibrating the simulated Cumulative Distribution Function (CDF)
#' by adding to the observed quantiles both the mean delta change and the individual delta changes in the corresponding quantiles.
#' This method is applicable to any kind of variable.
#'
#' It can keep the extreme value, if you choose constant extrapolation method. But then you will face the risk
#' that the extreme value is an error.
#'
#' \strong{gqm}
#'
#' Gamma Quantile Mapping. 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.
#'
#' It can somehow filter some extreme values caused by errors, while keep the extreme value. Seems more reasonable.
#' Better have a long period of training, and the if the forecast system is relatively stable.
#'
#' It is a generic function, if in your case you need to debug, please see \code{?debug()}
#' for how to debug S4 method.
#'
#' @examples
#'
#' ######## hyfo grid file biascorrection
#' ########
#'
#' # If your input is obtained by \code{loadNcdf}, you can also directly biascorrect
#' # the file.
#'
#' # First load ncdf file.
#' filePath <- system.file("extdata", "tnc.nc", package = "hyfo")
#' varname <- getNcdfVar(filePath)
#' nc <- loadNcdf(filePath, varname)
#'
#' data(tgridData)
#' # Since the example data, has some NA values, the process will include some warning #message,
#' # which can be ignored in this case.
#'
#'
#'
#'
#' # Then we will use nc data as forecasting data, and use itself as hindcast data,
#' # use tgridData as observation.
#' newFrc <- biasCorrect(nc, nc, tgridData)
#' newFrc <- biasCorrect(nc, nc, tgridData, scaleType = 'add')
#' newFrc <- biasCorrect(nc, nc, tgridData, method = 'eqm', extrapolate = 'constant',
#' preci = TRUE)
#' newFrc <- biasCorrect(nc, nc, tgridData, method = 'gqm', preci = TRUE)
#'
#'
#' ######## Time series biascorrection
#' ########
#'
#' # Use the time series from testdl as an example, we take frc, hindcast and obs from testdl.
#' data(testdl)
#'
#' # common period has to be extracted in order to better train the forecast.
#'
#' datalist <- extractPeriod(testdl, startDate = '1994-1-1', endDate = '1995-10-1')
#'
#' frc <- datalist[[1]]
#' hindcast <- datalist[[2]]
#' obs <- datalist[[3]]
#'
#'
#' # The data used here is just for example, so there could be negative data.
#'
#' # default method is scaling, with 'multi' scaleType
#' frc_new <- biasCorrect(frc, hindcast, obs)
#'
#' # for precipitation data, extra process needs to be executed, so you have to tell
#' # the program that it is a precipitation data.
#'
#' frc_new1 <- biasCorrect(frc, hindcast, obs, preci = TRUE)
#'
#' # You can use other scaling methods to biascorrect.
#' frc_new2 <- biasCorrect(frc, hindcast, obs, scaleType = 'add')
#'
#' #
#' frc_new3 <- biasCorrect(frc, hindcast, obs, method = 'eqm', preci = TRUE)
#' frc_new4 <- biasCorrect(frc, hindcast, obs, method = 'gqm', preci = TRUE)
#'
#' plotTS(obs, frc, frc_new, frc_new1, frc_new2, frc_new3, frc_new4, plot = 'cum')
#'
#' # You can also give name to this input list.
#' TSlist <- list(obs, frc, frc_new, frc_new1, frc_new2, frc_new3, frc_new4)
#' names(TSlist) <- c('obs', 'frc', 'delta', 'delta_preci', 'scale', 'eqm', 'gqm')
#' plotTS(list = TSlist, plot = 'cum')
#'
#'
#' # If the forecasts you extracted only has incontinuous data for certain months and years, e.g.,
#' # for seasonal forecasting, forecasts only provide 3-6 months data, so the case can be
#' # for example Dec, Jan and Feb of every year from year 1999-2005.
#' # In such case, you need to extract certain months and years from observed time series.
#' # extractPeriod() can be then used.
#'
#'
#'
#'
#'
#' # More examples can be found in the user manual on https://yuanchao-xu.github.io/hyfo/
#'
#'
#' @references
#' Bias correction methods come from \code{biasCorrection} from \code{dowscaleR}
#'
#' \itemize{
#'
#' \item Santander Meteorology Group (2015). downscaleR: Climate data manipulation and statistical downscaling. R
#' package version 0.6-0. https://github.com/SantanderMetGroup/downscaleR/wiki
#'
#' \item R.A.I. Wilcke, T. Mendlik and A. Gobiet (2013) Multi-variable error correction of regional climate models. Climatic Change, 120, 871-887
#'
#' \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 Yuanchao Xu \email{xuyuanchao37@@gmail.com }
#' @importFrom methods setMethod
#' @export
#'
setGeneric('biasCorrect', function(frc, hindcast, obs, method = 'scaling', scaleType = 'multi',
preci = FALSE, prThreshold = 0, extrapolate = 'no') {
standardGeneric('biasCorrect')
})
# Since in new version of roxygen2, describeIn was changed, https://stackoverflow.com/questions/24246594/automatically-document-all-methods-of-an-s4-generic-using-roxygen2
# so use rdname instead
#' @rdname biasCorrect
#'
setMethod('biasCorrect', signature('data.frame', 'data.frame', 'data.frame'),
function(frc, hindcast, obs, method, scaleType, preci, prThreshold, extrapolate) {
result <- biasCorrect.TS(frc, hindcast, obs, method, scaleType, preci, prThreshold, extrapolate)
return(result)
})
#' @rdname biasCorrect
setMethod('biasCorrect', signature('list', 'list', 'list'),
function(frc, hindcast, obs, method, scaleType, preci, prThreshold, extrapolate) {
result <- biasCorrect.list(frc, hindcast, obs, method, scaleType, preci, prThreshold, extrapolate)
return(result)
})
biasCorrect.TS <- function(frc, hindcast, obs, method, scaleType, preci, prThreshold, extrapolate) {
# First check if the first column is Date
if (!grepl('-|/', obs[1, 1]) | !grepl('-|/', hindcast[1, 1]) | !grepl('-|/', frc[1, 1])) {
stop('First column is not date or Wrong Date formate, check the format in ?as.Date{base}
and use as.Date to convert.If your input is a hyfo dataset, put input = "hyfo" as an
argument, check help for more info.')
}
# change to date type is easier, but in case in future the flood part is added, Date type doesn't have
# hour, min and sec, so, it's better to convert it into POSIxlt.
# if condition only accepts one condition, for list comparison, there are a lot of conditions, better
# further process it, like using any.
if (any(as.POSIXlt(hindcast[, 1]) != as.POSIXlt(obs[, 1]))) {
warning('time of obs and time of hindcast are not the same, which may cause inaccuracy in
the calibration.')
}
n <- ncol(frc)
# For every column, it's biascorrected respectively.
frc_data <- lapply(2:n, function(x) biasCorrect_core(frc[, x], hindcast[, x], obs[, 2], method = method,
scaleType = scaleType, preci = preci, prThreshold = prThreshold,
extrapolate = extrapolate))
frc_data <- do.call('cbind', frc_data)
rownames(frc_data) <- NULL
names <- colnames(frc)
frc_new <- data.frame(frc[, 1], frc_data)
colnames(frc_new) <- names
return(frc_new)
}
biasCorrect.list <- function(frc, hindcast, obs, method, scaleType, preci, prThreshold, extrapolate) {
## Check if the data is a hyfo grid data.
checkHyfo(frc, hindcast, obs)
hindcastData <- hindcast$Data
obsData <- obs$Data
frcData <- frc$Data
## save frc dimension order, at last, set the dimension back to original dimension
frcDim <- attributes(frcData)$dimensions
## ajust the dimension into general dimension order.
hindcastData <- adjustDim(hindcastData, ref = c('lon', 'lat', 'time'))
obsData <- adjustDim(obsData, ref = c('lon', 'lat', 'time'))
## CheckDimLength, check if all the input dataset has different dimension length
# i.e. if they all have the same lon and lat number.
checkDimLength(frcData, hindcastData, obsData, dim = c('lon', 'lat'))
# Now real bias correction is executed.
# for some forcasts, they have results from different models or scenarios, if so
# there will be a dimension called member
memberIndex <- grepAndMatch('member', attributes(frcData)$dimensions)
# For dataset that has a member part
if (length(memberIndex) != 0) {
# check if frcData and hindcastData has the same dimension and length.
checkDimLength(frcData, hindcastData, dim = 'member')
frcData <- adjustDim(frcData, ref = c('lon', 'lat', 'time', 'member'))
# The following code may speed up because it doesn't use for loop.
# It firstly combine different array into one array. combine the time
# dimension of frc, hindcast and obs. Then use apply, each time extract
# the total time dimension, and first part is frc, second is hindcast, third
# is obs. Then use these three parts to bias correct. All above can be written
# in one function and called within apply. But too complicated to understand,
# So save it for future use maybe.
# for (member in 1:dim(frcData)[4]) {
# totalArr <- array(c(frcData[,,, member], hindcastData[,,, member], obsData),
# dim = c(dim(frcData)[1], dim(frcData)[2],
# dim(frcData)[3] + dim(hindcastData)[3] + dim(obsData)[3]))
# }
for (member in 1:dim(frcData)[4]) {
for (lon in 1:dim(frcData)[1]) {
for (lat in 1:dim(frcData)[2]) {
frcData[lon, lat,, member] <- biasCorrect_core(frcData[lon, lat,,member], hindcastData[lon, lat,, member], obsData[lon, lat,], method = method,
scaleType = scaleType, preci = preci, prThreshold = prThreshold,
extrapolate = extrapolate)
}
}
}
} else {
frcData <- adjustDim(frcData, ref = c('lon', 'lat', 'time'))
for (lon in 1:dim(frcData)[1]) {
for (lat in 1:dim(frcData)[2]) {
frcData[lon, lat,] <- biasCorrect_core(frcData[lon, lat,], hindcastData[lon, lat,], obsData[lon, lat,], method = method,
scaleType = scaleType, preci = preci, prThreshold = prThreshold,
extrapolate = extrapolate)
}
}
}
frcData <- adjustDim(frcData, ref = frcDim)
frc$Data <- frcData
frc$biasCorrected_by <- method
frc_new <- frc
return(frc_new)
}
#' @importFrom MASS fitdistr
#' @importFrom stats ecdf quantile pgamma qgamma rgamma
#'
#' @references
#' Bias correction methods come from \code{biasCorrection} from \code{dowscaleR}
#'
#' \itemize{
#'
#' \item Santander Meteorology Group (2015). downscaleR: Climate data manipulation and statistical downscaling. R
#' package version 0.6-0. https://github.com/SantanderMetGroup/downscaleR/wiki
#'
#' \item R.A.I. Wilcke, T. Mendlik and A. Gobiet (2013) Multi-variable error correction of regional climate models. Climatic Change, 120, 871-887
#'
#' \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
#' }
#'
#'
#'
# this is only used to calculate the value column,
biasCorrect_core <- function(frc, hindcast, obs, method, scaleType, preci, prThreshold, extrapolate){
# If the variable is precipitation, some further process needs to be added.
# The process is taken from downscaleR, to provide a more reasonable hindcast, used in the calibration.
# check if frc, hindcast or obs are all na values
if (!any(!is.na(obs)) | !any(!is.na(frc)) | !any(!is.na(hindcast))) {
warning('In this cell, frc, hindcast or obs data is missing. No biasCorrection for this cell.')
return(NA)
}
if (preci == TRUE) {
preprocessHindcast_res <- preprocessHindcast(hindcast = hindcast, obs = obs, prThreshold = prThreshold)
hindcast <- preprocessHindcast_res[[1]]
minHindcastPreci <- preprocessHindcast_res[[2]]
}
# default is the simplest method in biascorrection, just do simple addition and subtraction.
if (method == 'delta') {
if (length(frc) != length(obs)) stop('This method needs frc data have the same length as obs data.')
# comes from downscaleR biascorrection method
frcMean <- mean(frc, na.rm = TRUE)
hindcastMean <- mean(hindcast, na.rm = TRUE)
frc <- obs - hindcastMean + frcMean
} else if (method == 'scaling') {
obsMean <- mean(obs, na.rm = TRUE)
hindcastMean <- mean(hindcast, na.rm = TRUE)
if (scaleType == 'multi') {
frc <- frc / hindcastMean * obsMean
} else if (scaleType == 'add') {
frc <- frc - hindcastMean + obsMean
}
} else if (method == 'eqm') {
if (preci == FALSE) {
frc <- biasCorrect_core_eqm_nonPreci(frc, hindcast, obs, extrapolate, prThreshold)
} else {
frc <- biasCorrect_core_eqm_preci(frc, hindcast, obs, minHindcastPreci, extrapolate,
prThreshold)
}
} else if (method == 'gqm') {
if (preci == FALSE) stop ('gqm method only applys to precipitation, please set preci = T')
frc <- biasCorrect_core_gqm(frc, hindcast, obs, prThreshold, minHindcastPreci)
}
return(frc)
}
#' @importFrom MASS fitdistr
#' @importFrom stats rgamma
preprocessHindcast <- function(hindcast, obs, prThreshold) {
lowerIndex <- length(which(obs < prThreshold))
# In the original function, this minHindcastPreci is Pth[,i,j] in downscaleR, and it is originally
# set to NA, which is not so appropriate for all the precipitations.
# In the original function, there are only two conditions, 1. all the obs less than threshold
# 2. there are some obs less than threshold.
# While, if we set threshold to 0, there could be a 3rd condition, all the obs no less than threshold.
# Here I set this situation, firstly set minHindcastPreci to the min of the hindcast. Because in future
# use, 'eqm' method is going to use this value.
# The problem above has been solved.
if (lowerIndex >= 0 & lowerIndex < length(obs)) {
index <- sort(hindcast, decreasing = FALSE, na.last = NA, index.return = TRUE)$ix
hindcast_sorted <- sort(hindcast, decreasing = FALSE, na.last = NA)
# minHindcastPreci is the min preci over threshold FOR ***HINDCAST***
# But use obs to get the lowerIndex, so obs_sorted[lowerIndex + 1] > prThreshold, but
# hindcast_sorted[lowerIndex + 1] may greater than or smaller than ptThreshold
# It would be better to understand if you draw two lines: hindcast_sorted and obs_sorted
# with y = prThreshold, you will find the difference of the two.
# In principle, the value under the threshold needs to be replaced by some other reasonable value.
# simplest way
minHindcastPreci <- hindcast_sorted[lowerIndex + 1]
# Also here if minHindcastPreci is 0 and prThreshold is 0, will cause problem, bettter set
# I set it prThreshold != 0
if (minHindcastPreci <= prThreshold & prThreshold != 0) {
obs_sorted <- sort(obs, decreasing = FALSE, na.last = NA)
# higherIndex is based on hindcast
higherIndex <- which(hindcast_sorted > prThreshold & !is.na(hindcast_sorted))
if (length(higherIndex) == 0) {
higherIndex <- max(which(!is.na(hindcast_sorted)))
higherIndex <- min(length(obs_sorted), higherIndex)
} else {
higherIndex <- min(higherIndex)
}
# here I don't know why choose 6.
# Written # [Shape parameter Scale parameter] in original package
# according to the reference and gamma distribution, at least 6 values needed to fit gamma
# distribution.
if (length(unique(obs_sorted[(lowerIndex + 1):higherIndex])) < 6) {
hindcast_sorted[(lowerIndex + 1):higherIndex] <- mean(obs_sorted[(lowerIndex + 1):higherIndex],
na.rm = TRUE)
} else {
obsGamma <- fitdistr(obs_sorted[(lowerIndex + 1):higherIndex], "gamma")
# this is to replace the original hindcast value between lowerIndex and higherIndex with
# some value taken from gamma distribution just generated.
hindcast_sorted[(lowerIndex + 1):higherIndex] <- rgamma(higherIndex - lowerIndex, obsGamma$estimate[1],
rate = obsGamma$estimate[2])
}
hindcast_sorted <- sort(hindcast_sorted, decreasing = FALSE, na.last = NA)
}
minIndex <- min(lowerIndex, length(hindcast))
hindcast_sorted[1:minIndex] <- 0
hindcast[index] <- hindcast_sorted
} else if (lowerIndex == length(obs)) {
index <- sort(hindcast, decreasing = FALSE, na.last = NA, index.return = TRUE)$ix
hindcast_sorted <- sort(hindcast, decreasing = FALSE, na.last = NA)
minHindcastPreci <- hindcast_sorted[lowerIndex]
# here is to compare with hindcast, not obs
minIndex <- min(lowerIndex, length(hindcast))
hindcast_sorted[1:minIndex] <- 0
hindcast[index] <- hindcast_sorted
}
return(list(hindcast, minHindcastPreci))
}
biasCorrect_core_eqm_nonPreci <- function(frc, hindcast, obs, extrapolate, prThreshold) {
ecdfHindcast <- ecdf(hindcast)
if (extrapolate == 'constant') {
higherIndex <- which(frc > max(hindcast, na.rm = TRUE))
lowerIndex <- which(frc < min(hindcast, na.rm = TRUE))
extrapolateIndex <- c(higherIndex, lowerIndex)
non_extrapolateIndex <- setdiff(1:length(frc), extrapolateIndex)
# In case extrapolateIndex is of length zero, than extrapolate cannot be used afterwards
# So use setdiff(1:length(sim), extrapolateIndex), if extrapolateIndex == 0, than it will
# return 1:length(sim)
if (length(higherIndex) > 0) {
maxHindcast <- max(hindcast, na.rm = TRUE)
dif <- maxHindcast - max(obs, na.rm = TRUE)
frc[higherIndex] <- frc[higherIndex] - dif
}
if (length(lowerIndex) > 0) {
minHindcast <- min(hindcast, na.rm = TRUE)
dif <- minHindcast - min(obs, nna.rm = TRUE)
frc[lowerIndex] <- frc[lowerIndex] - dif
}
frc[non_extrapolateIndex] <- quantile(obs, probs = ecdfHindcast(frc[non_extrapolateIndex]),
na.rm = TRUE, type = 4)
} else {
frc <- quantile(obs, probs = ecdfHindcast(frc), na.rm = TRUE, type = 4)
}
return(frc)
}
biasCorrect_core_eqm_preci <- function(frc, hindcast, obs, minHindcastPreci, extrapolate,
prThreshold) {
# Most of time this condition seems useless because minHindcastPreci comes from hindcast, so there will be
# always hindcast > minHindcastPreci exists.
# Unless one condition that minHindcastPreci is the max in the hindcast, than on hindcast > minHindcastPreci
if (length(which(hindcast > minHindcastPreci)) > 0) {
ecdfHindcast <- ecdf(hindcast[hindcast > minHindcastPreci])
noRain <- which(frc <= minHindcastPreci & !is.na(frc))
rain <- which(frc > minHindcastPreci & !is.na(frc))
# drizzle is to see whether there are some precipitation between the min frc (over threshold) and
# min hindcast (over threshold).
drizzle <- which(frc > minHindcastPreci & frc <= min(hindcast[hindcast > minHindcastPreci], na.rm = TRUE)
& !is.na(frc))
if (length(rain) > 0) {
ecdfFrc <- ecdf(frc[rain])
if (extrapolate == 'constant') {
# This higher and lower index mean the extrapolation part
higherIndex <- which(frc[rain] > max(hindcast, na.rm = TRUE))
lowerIndex <- which(frc[rain] < min(hindcast, na.rm = TRUE))
extrapolateIndex <- c(higherIndex, lowerIndex)
non_extrapolateIndex <- setdiff(1:length(rain), extrapolateIndex)
if (length(higherIndex) > 0) {
maxHindcast <- max(hindcast, na.rm = TRUE)
dif <- maxHindcast - max(obs, na.rm = TRUE)
frc[rain[higherIndex]] <- frc[higherIndex] - dif
}
if (length(lowerIndex) > 0) {
minHindcast <- min(hindcast, na.rm = TRUE)
dif <- minHindcast - min(obs, nna.rm = TRUE)
frc[rain[lowerIndex]] <- frc[lowerIndex] - dif
}
# Here the original function doesn't accout for the situation that extraploateIndex is 0
# if it is 0, rain[-extraploateIndex] would be nothing
# Above has been solved by using setdiff.
frc[rain[non_extrapolateIndex]] <- quantile(obs[which(obs > prThreshold & !is.na(obs))],
probs = ecdfHindcast(frc[rain[non_extrapolateIndex]]),
na.rm = TRUE, type = 4)
} else {
frc[rain] <- quantile(obs[which(obs > prThreshold & !is.na(obs))],
probs = ecdfHindcast(frc[rain]), na.rm = TRUE, type = 4)
}
}
if (length(drizzle) > 0){
# drizzle part is a seperate part. it use the ecdf of frc (larger than minHindcastPreci) to
# biascorrect the original drizzle part
frc[drizzle] <- quantile(frc[which(frc > min(hindcast[which(hindcast > minHindcastPreci)], na.rm = TRUE) &
!is.na(frc))], probs = ecdfFrc(frc[drizzle]), na.rm = TRUE,
type = 4)
}
frc[noRain] <- 0
} else {
# in this condition minHindcastPreci is the max of hindcast, so all hindcast <= minHindcastPreci
# And frc distribution is used then.
noRain <- which(frc <= minHindcastPreci & !is.na(frc))
rain <- which(frc > minHindcastPreci & !is.na(frc))
if (length(rain) > 0) {
ecdfFrc <- ecdf(frc[rain])
frc[rain] <- quantile(obs[which(obs > prThreshold & !is.na(obs))], probs = ecdfFrc(frc[rain]),
na.rm = TRUE, type = 4)
}
frc[noRain]<-0
}
return(frc)
}
biasCorrect_core_gqm <- function(frc, hindcast, obs, prThreshold, minHindcastPreci) {
if (any(obs > prThreshold)) {
ind <- which(obs > prThreshold & !is.na(obs))
obsGamma <- fitdistr(obs[ind],"gamma")
ind <- which(hindcast > 0 & !is.na(hindcast))
hindcastGamma <- fitdistr(hindcast[ind],"gamma")
rain <- which(frc > minHindcastPreci & !is.na(frc))
noRain <- which(frc <= minHindcastPreci & !is.na(frc))
probF <- pgamma(frc[rain], hindcastGamma$estimate[1], rate = hindcastGamma$estimate[2])
frc[rain] <- qgamma(probF,obsGamma$estimate[1], rate = obsGamma$estimate[2])
frc[noRain] <- 0
} else {
warning('All the observations of this cell(station) are lower than the threshold,
no bias correction applied.')
}
return(frc)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.