2A6E2BEA

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