diff --git a/.Rhistory b/.Rhistory index 526c63d..5ed785c 100644 --- a/.Rhistory +++ b/.Rhistory @@ -1,325 +1,12 @@ } -if (cell == 'mean') { -TS <- apply(data, MARGIN = 3, FUN = mean, na.rm = TRUE) -} else { -TS <- data[cell[1], cell[2], ] -} -return(TS) -} -#' @importFrom reshape2 melt -#' @describeIn getPreciBar -getPreciBar.TS <- function(TS) { -Date <- as.POSIXlt(TS[, 1]) -yearIndex <- Date$year + 1900 -monthIndex <- Date$mon + 1 -n <- ncol(TS) - 1 -if ( n == 1) { -TS <- TS[, 2] -} else { -TS <- TS[, -1] -# month index should be repeat, but years cannot. -yearIndex <- sapply(1:n, function(x) yearIndex + x - 1) -dim(yearIndex) <- c(n * nrow(yearIndex), 1) -monthIndex <- rep(monthIndex, n) -TS <- melt(TS)[, 2] -} -return(TS) -} -#' @importFrom stats median -#' @importFrom reshape2 melt -#' @import ggplot2 -#' @describeIn getPreciBar -getPreciBar.plot <- function(TS, method, output, name, plotRange, omitNA, info, -yearIndex, monthIndex, ...) { -if (method == 'meanMonthly') { -monthlyPreci <- tapply(TS, INDEX = list(yearIndex, monthIndex), FUN = sum, na.rm = omitNA) -meanMonthlyPreci <- apply(monthlyPreci, MARGIN = 2, FUN = mean, na.rm = TRUE) -title <- 'Mean Monthly Precipitation' -xlab <- 'Month' -plotPreci <- data.frame(Index = month.abb[as.numeric(colnames(monthlyPreci))], -Preci = meanMonthlyPreci) -# Here factor has to be reassigned, to keep the original order, or it will be reordered. -plotPreci$Index <- factor(plotPreci$Index, levels = plotPreci$Index, ordered = TRUE) -if (plotRange) { -maxValue <- apply(monthlyPreci, MARGIN = 2, FUN = max, na.rm = TRUE) -minValue <- apply(monthlyPreci, MARGIN = 2, FUN = min, na.rm = TRUE) -plotPreci$maxValue <- maxValue -plotPreci$minValue <- minValue -ylim <- c(0,max(maxValue, na.rm = TRUE) * 1.1) -} else { -ylim <- c(0,max(meanMonthlyPreci, na.rm = TRUE) * 1.1) -} -} else if (method == 'annual') { -if (length(unique(monthIndex)) < 12) { -warning ('There are less than 12 months in a year, the results may be inaccurate.') -} -annualPreci <- tapply(TS, INDEX = yearIndex, FUN = sum, na.rm = TRUE) -title <- 'Annual Precipitation' -xlab <- 'Year' -plotName <- names(annualPreci) -plotPreci <- data.frame(Index = names(annualPreci), Preci = annualPreci) -plotPreci$Index <- factor(plotPreci$Index, levels = plotPreci$Index, ordered = TRUE) -ylim <- c(0, max(annualPreci, na.rm = TRUE) * 1.1) -} else if (is.numeric(method)) { -month <- method -monExisting <- length(which(unique(monthIndex) == month)) -if (monExisting == 0) stop("Your input month doesn't exist in the dataset.") -monthlyPreci <- getMeanPreci(TS, method = month, yearIndex = yearIndex, -monthIndex = monthIndex, fullResults = TRUE, omitNA = omitNA) -# If monthlyPreci length is 1, names need to be added. -if (length(monthlyPreci) == 1) names(monthlyPreci) <- unique(yearIndex) -plotPreci <- data.frame(Index = names(monthlyPreci), Preci = monthlyPreci) -plotPreci$Index <- factor(plotPreci$Index, levels = plotPreci$Index, ordered = TRUE) -title <- paste(month.abb[month], 'Precipitation over Whole Period', sep = ' ') -xlab <- 'Year' -ylim <- c(0, max(monthlyPreci, na.rm = TRUE) * 1.1) -} else if (method == 'spring') { -wm <- match(c(3, 4, 5), unique(monthIndex)) -if (length(which(!is.na(wm))) < 3) { -stop('Spring has less than 3 months, check data and try to calculate every month -seperately or choose another season.') -} -seasonalPreci <- getMeanPreci(TS, method = 'spring', yearIndex = yearIndex, -monthIndex = monthIndex, fullResults = TRUE, omitNA = omitNA) -plotPreci <- data.frame(Index = names(seasonalPreci), Preci = seasonalPreci) -plotPreci$Index <- factor(plotPreci$Index, levels = plotPreci$Index, ordered = TRUE) -title <- paste('Spring', 'Precipitation over Whole Period', sep = ' ') -xlab <- 'Year' -ylim <- c(0, max(seasonalPreci, na.rm = TRUE) * 1.1) -} else if (method == 'summer') { -wm <- match(c(6, 7, 8), unique(monthIndex)) -if (length(which(!is.na(wm))) < 3) { -stop('Summer has less than 3 months, check data and try to calculate every month -seperately or choose another season.') -} -seasonalPreci <- getMeanPreci(TS, method = 'summer', yearIndex = yearIndex, -monthIndex = monthIndex, fullResults = TRUE, omitNA = omitNA) -plotPreci <- data.frame(Index = names(seasonalPreci), Preci = seasonalPreci) -plotPreci$Index <- factor(plotPreci$Index, levels = plotPreci$Index, ordered = TRUE) -title <- paste('Summer', 'Precipitation over Whole Period', sep = ' ') -xlab <- 'Year' -ylim <- c(0, max(seasonalPreci, na.rm = TRUE) * 1.1) -} else if (method == 'autumn') { -wm <- match(c(9, 10, 11), unique(monthIndex)) -if (length(which(!is.na(wm))) < 3) { -stop('Autumn has less than 3 months, check data and try to calculate every month -seperately or choose another season.') -} -seasonalPreci <- getMeanPreci(TS, method = 'autumn', yearIndex = yearIndex, -monthIndex = monthIndex, fullResults = TRUE, omitNA = omitNA) -plotPreci <- data.frame(Index = names(seasonalPreci), Preci = seasonalPreci) -plotPreci$Index <- factor(plotPreci$Index, levels = plotPreci$Index, ordered = TRUE) -title <- paste('Autumn', 'Precipitation over Whole Period', sep = ' ') -xlab <- 'Year' -ylim <- c(0, max(seasonalPreci, na.rm = TRUE) * 1.1) -} else if (method == 'winter') { -wm <- match(c(12, 1, 2), unique(monthIndex)) -if (length(which(!is.na(wm))) < 3) { -stop('Winter has less than 3 months, check data and try to calculate every month -seperately or choose another season.') -} -seasonalPreci <- getMeanPreci(TS, method = 'winter', yearIndex = yearIndex, -monthIndex = monthIndex, fullResults = TRUE, omitNA = omitNA) -plotPreci <- data.frame(Index = names(seasonalPreci), Preci = seasonalPreci) -plotPreci$Index <- factor(plotPreci$Index, levels = plotPreci$Index, ordered = TRUE) -title <- paste('Winter', 'Precipitation over Whole Period', sep = ' ') -xlab <- 'Year' -ylim <- c(0, max(seasonalPreci, na.rm = TRUE) * 1.1) -} else { -stop(paste('No method called "', method, '", check help for information')) -} -xlim <- c(0, length(rownames(plotPreci))) -if (info == TRUE) { -meanValue <- round(mean(plotPreci$Preci, na.rm = TRUE), 2) -medianValue <- round(median(plotPreci$Preci,na.rm = TRUE), 2) -plotMean <- paste('Mean', ' = ', meanValue) -plotMedian <- paste('Median', ' = ', medianValue) -plotMax <- round(max(plotPreci$Preci, na.rm = TRUE), 2) -plotMin <- round(min(plotPreci$Preci, na.rm = TRUE), 2) -word <- paste('\n\n', paste(' Max', '=', plotMax), ',', paste('Min', '=', plotMin), ',', -plotMean, ',', plotMedian) -} else word <- NULL -xlab <- paste(xlab, word) -theme_set(theme_bw()) -mainLayer <- with(plotPreci, { -ggplot(plotPreci) + -geom_bar(aes(x = Index, y = Preci), stat = 'identity', colour = 'black', fill = 'cyan2', width = rel(.4)) + -xlab(xlab) + -ylab('Precipitation (mm)') + -ggtitle(title) + -labs(empty = NULL, ...) +#in order to pass "...", arguments shouldn't be empty. -theme(plot.title = element_text(size = rel(1.6), face = 'bold'), -axis.title.x = element_text(size = rel(1.6)), -axis.title.y = element_text(size = rel(1.6)), -axis.text.x = element_text(angle = 90, hjust = 1, size = rel(1.9)), -axis.text.y = element_text(size = rel(1.9))) -# geom_text(x = min(xlim) + 0.95 * (max(xlim) - min(xlim)), y = min(ylim) + 0.15 * (max(ylim) - min(ylim)), -# label = word)+ -# geom_hline(yintercept = meanValue) + -# geom_text(x = min(xlim) + 0.3 * (max(xlim) - min(xlim)), y = meanValue + 3, vjust = 0, label = 'mean') + -# geom_hline(yintercept = medianValue, colour = 'red') + -# geom_text(x = min(xlim) + 0.6 * (max(xlim) - min(xlim)), y = medianValue + 3, vjust = 0, -# label = 'median', colour = 'red') -}) -if (plotRange) { -if (is.null(plotPreci$maxValue)) { -message('There is no plotRange for this method') -print(mainLayer) -} else { -rangeLayer <- with(plotPreci, { -geom_errorbar(aes(x = Index, ymax = maxValue, ymin = minValue), width = rel(0.3)) -}) -print(mainLayer + rangeLayer) -} -} else { -print(mainLayer) -} -if (output == 'plot') { -return(mainLayer) -} else if (output == 'ggplot') { -if (is.null(name)) stop('"name" argument not found, -If you choose "ggplot" as output, please assign a name.') -plotPreci$Name <- rep(name, dim(plotPreci)[1]) -return(plotPreci) -} else { -return(plotPreci) -} -} -#' Combine bars together -#' @param ... different barplots generated by \code{getPreciBar(, output = 'ggplot')}, refer to details. -#' @details -#' ..., representing different ouput generated by \code{getPreciBar(, output = 'ggplot')}, they -#' have to be of the same type, e.g., -#' 1. Jan precipitation of different years, Feb precipitation of different years, and... -#' They are both monthly precipitation, and they share x axis. -#' -#' 2. Mean monthly precipitation of different dataset. e.g., long term mean monthly precipitation -#' and short term mean monthly precipitation. They are both mean monthly precipitation. -#' -#' @param nrow A number showing the number of rows. -#' @param list If input is a list containing different ggplot data, use l\code{list = inputlist}. -#' NOTE: yOU HAVE TO PUT A \code{list = }, before your list. -#' @param x A string of x axis name. -#' @param y A string of y axis name. -#' @param title A string of the title. -#' @param output A boolean, if chosen TRUE, the output will be given. -#' @return A combined barplot. -#' @examples -#' -#' data(tgridData)# the result of \code{loadGridData{ecomsUDG.Raccess}} -#' #output type of getPreciBar() has to be 'ggplot'. -#' b1 <- getPreciBar(tgridData, method = 2, output = 'ggplot', name = 'b1') -#' b2 <- getPreciBar(tgridData, method = 3, output = 'ggplot', name = 'b2') -#' -#' getPreciBar_comb(b1, b2) -#' -#' # More examples can be found in the user manual on http://yuanchao-xu.github.io/hyfo/ -#' -#' @export -#' @import ggplot2 -#' @references -#' -#' \itemize{ -#' \item H. Wickham. ggplot2: elegant graphics for data analysis. Springer New York, 2009. -#' } -#' -getPreciBar_comb <- function(..., list = NULL, nrow = 1, x = '', y = '', title = '', output = FALSE) { -if (!is.null(list)) { -data_ggplot <- do.call('rbind', list) -} else { -bars <- list(...) -checkBind(bars, 'rbind') -data_ggplot <- do.call('rbind', bars) -} -if (!class(data_ggplot) == 'data.frame') { -warning('Your input is probably a list, but you forget to add "list = " before it. -Try again, or check help for more information.') -} else if (is.null(data_ggplot$Name)) { -stop('No "Name" column in the input data, check the arguments in getPreciBar(), if -output = "ggplot" is assigned, more info please check ?getPreciBar.') -} -data_ggplot$Name <- factor(data_ggplot$Name, levels = unique(data_ggplot$Name), ordered = TRUE) -theme_set(theme_bw()) -mainLayer <- with(data_ggplot, { -ggplot(data_ggplot) + -geom_bar(aes(x = Index, y = Preci),fill = 'cyan2', stat = 'identity', -colour = 'black', width = rel(.4)) + -facet_wrap( ~ Name, nrow = nrow) + -theme(plot.title = element_text(size = rel(1.6), face = 'bold'), -axis.title.x = element_text(size = rel(1.6)), -axis.title.y = element_text(size = rel(1.6)), -axis.text.x = element_text(angle = 90, hjust = 1, size = rel(1.9)), -axis.text.y = element_text(size = rel(1.9))) + -labs(x = x, y = y, title = title) -}) -if (!any(is.na(match(c('minValue', 'maxValue'), colnames(data_ggplot))))) { -rangeLayer <- with(data_ggplot, { -geom_errorbar(aes(x = Index, ymax = maxValue, ymin = minValue), width = rel(0.3)) -}) -mainLayer <- mainLayer + rangeLayer -} -suppressWarnings(print(mainLayer)) -if (output == TRUE) return(data_ggplot) -} -b1 <- getPreciBar(tgridData, method = 'annual') -devtools::document() -devtools::document() -devtools::document() -devtools::check() -devtools::document() -devtools::check() -devtools::document() -devtools::check() -b1 <- getPreciBar(tgridData, method = 'annual') -devtools::check() -b1 <- getPreciBar(tgridData, method = 'annual') -b1 <- getPreciBar(tgridData, method = 'annual', cell(4,5)) -b1 <- getPreciBar(tgridData, method = 'annual', cell = c(4,5)) -Q -devtools::check() -b1 <- getPreciBar(tgridData, method = 'annual', cell = c(4,5)) -b1 <- getPreciBar(tgridData, method = 'annual', cell = c(14,5)) -b1 <- getPreciBar(tgridData, method = 'annual', cell = c(10,5)) -b1 <- getPreciBar(tgridData, method = 'annual', cell = c(10,2)) -b1 <- getPreciBar(tgridData, method = 'annual', cell = c(5,5)) -debug(getPreciBar) -b1 <- getPreciBar(tgridData, method = 'annual', cell = c(5,5)) -n -Q -debug(getPreciBar.list) -b1 <- getPreciBar(tgridData, method = 'annual', cell = c(5,5)) -undebug(getPreciBar) -b1 <- getPreciBar(tgridData, method = 'annual', cell = c(5,5)) -data[cell[1], cell[2], ] -str(data) -data[5, 3, ] -b1 <- getPreciBar(tgridData, method = 'annual', cell = c(5,3)) -b1 <- getPreciBar(tgridData, method = 'annual') -b1 <- getPreciBar(tgridData, method = 'annual', cell = c(5,3)) -str(data) -b1 <- getPreciBar(tgridData, method = 2, output = 'ggplot', name = 'b1') -undebug(getPreciBar.list) -b2 <- getPreciBar(tgridData, method = 3, output = 'ggplot', name = 'b2') -b1 <- getPreciBar(tgridData, method = 2, output = 'ggplot', name = 'b1') -getPreciBar_comb(b1, b2) -getPreciBar_comb(b1, b2, nrow = 2) -devtools::document() -devtools::document() -devtools::check() -devtools::document() -devtools::document() -devtools::document() -devtools::document() -devtools::document() -devtools::document() -devtools::document() -devtools::document() -devtools::document() -devtools::document() -devtools::check() -devtools::document() -devtools::check() -debug(hyfoUpdates) +.onAttach <- function(libname, pkgname) { +message_out <- suppressWarnings(try(hyfoUpdates(), silent = TRUE)) +if (!is.null(message_out)) { +if (grepl('Version', message_out)) { +packageStartupMessage(message_out) +} +} +} ## For package updates information #' @importFrom utils packageDescription hyfoUpdates <- function(){ @@ -328,13 +15,17 @@ updatesLine <- grep('id=\\"updates"', page) versionLine <- updatesLine + 2 version <- unlist(strsplit(page[versionLine], split = ' '))[2] version_local <- packageDescription("hyfo")$Version +# the first tow digit is the most important part of the version +version12 <- unlist(strsplit(version, split = "[.]"))[1:2] +version_local12 <- unlist(strsplit(version_local1, split = "[.]"))[1:2] +sameVersion <- version12 == version_local12 # generate message version_msg <- strsplit(strsplit(page[versionLine], split = '

')[[1]][2], split = '

')[[1]] infoLine <- versionLine + 2 info_msg <- strsplit(strsplit(page[infoLine], split = '

')[[1]][2], split = '

')[[1]] install_msg <- 'You can update by type in: devtools::install_gihub("Yuanchao-Xu/hyfo")' message_out <- NULL -if (version != version_local) { +if (!sameVersion) { message_out <- paste(version_msg, info_msg, install_msg, sep = '\n') } return(message_out) @@ -349,27 +40,6 @@ packageStartupMessage(message_out) } debug(hyfoUpdates) hyfoUpdates() -bn -page -upn -version -version_local -version_msg -strsplit('1.2.8', sep = '.') -strsplit('1.2.8') -strsplit('1.2.8', split = '.') -strsplit('1.2.8', split = '.')[1] -strsplit('1.2.8', split = '2') -strsplit('1.2.8', split = '.') -strsplit('1.2.8', split = '/.') -strsplit('1.2.8', split = '\.') -strsplit('1.2.8', split = '//.') -strsplit('1.2.8', split = "") -strsplit('1.2.8', split = ".") -strsplit('1.2.8', split = " .") -?strsplit -strsplit('1.2.8', split = "[.]") -unlist(strsplit('1.2.8', split = "[.]")) ## For package updates information #' @importFrom utils packageDescription hyfoUpdates <- function(){ @@ -380,7 +50,7 @@ version <- unlist(strsplit(page[versionLine], split = ' '))[2] version_local <- packageDescription("hyfo")$Version # the first tow digit is the most important part of the version version12 <- unlist(strsplit(version, split = "[.]"))[1:2] -version_local12 <- unlist(strsplit(version_local1, split = "[.]"))[1:2] +version_local12 <- unlist(strsplit(version_local, split = "[.]"))[1:2] sameVersion <- version12 == version_local12 # generate message version_msg <- strsplit(strsplit(page[versionLine], split = '

')[[1]][2], split = '

')[[1]] @@ -401,6 +71,11 @@ packageStartupMessage(message_out) } } } +debug(hyfoUpdates) +hyfoUpdates() +version12 +version_local12 +sameVersion ## For package updates information #' @importFrom utils packageDescription hyfoUpdates <- function(){ @@ -411,7 +86,7 @@ version <- unlist(strsplit(page[versionLine], split = ' '))[2] version_local <- packageDescription("hyfo")$Version # the first tow digit is the most important part of the version version12 <- unlist(strsplit(version, split = "[.]"))[1:2] -version_local12 <- unlist(strsplit(version_local1, split = "[.]"))[1:2] +version_local12 <- unlist(strsplit(version_local, split = "[.]"))[1:2] sameVersion <- version12 == version_local12 # generate message version_msg <- strsplit(strsplit(page[versionLine], split = '

')[[1]][2], split = '

')[[1]] @@ -419,7 +94,7 @@ infoLine <- versionLine + 2 info_msg <- strsplit(strsplit(page[infoLine], split = '

')[[1]][2], split = '

')[[1]] install_msg <- 'You can update by type in: devtools::install_gihub("Yuanchao-Xu/hyfo")' message_out <- NULL -if (!sameVersion) { +if (any(sameVersion == FALSE)) { message_out <- paste(version_msg, info_msg, install_msg, sep = '\n') } return(message_out) @@ -434,6 +109,20 @@ packageStartupMessage(message_out) } debug(hyfoUpdates) hyfoUpdates() +N +sameVersion +message_out +devtools::document() +devtools::check() +devtools::check() +library(hyfo) +devtools::install_github('Yuanchao-Xu/hyfo') +library(hyfo) +?getBiasFactor +?applyBiasFactor +?resample +library(hyfo) +library(hyfo) ## For package updates information #' @importFrom utils packageDescription hyfoUpdates <- function(){ @@ -452,7 +141,7 @@ infoLine <- versionLine + 2 info_msg <- strsplit(strsplit(page[infoLine], split = '

')[[1]][2], split = '

')[[1]] install_msg <- 'You can update by type in: devtools::install_gihub("Yuanchao-Xu/hyfo")' message_out <- NULL -if (!sameVersion) { +if (any(sameVersion == FALSE)) { message_out <- paste(version_msg, info_msg, install_msg, sep = '\n') } return(message_out) @@ -465,11 +154,80 @@ packageStartupMessage(message_out) } } } +hyfoUpdates() +debug(hyfoUpdates) +hyfoUpdates() +version_msg +info_msg +info_msg +strsplit(strsplit(page[infoLine], split = '

')[[1]][2], split = '

')[[1]] +strsplit(strsplit(page[infoLine + 1], split = '

')[[1]][2], split = '

')[[1]] +infoLine +page[57] +page[58] +page[58 : 60] +page[57] +page[60] +page[61] +page[63] +page[62] +page[64] +## For package updates information +#' @importFrom utils packageDescription +hyfoUpdates <- function(){ +page <- readLines('http://yuanchao-xu.github.io/hyfo/') +updatesLine <- grep('id=\\"updates"', page) +versionLine <- updatesLine + 2 +version <- unlist(strsplit(page[versionLine], split = ' '))[2] +version_local <- packageDescription("hyfo")$Version +# the first tow digit is the most important part of the version +version12 <- unlist(strsplit(version, split = "[.]"))[1:2] +version_local12 <- unlist(strsplit(version_local, split = "[.]"))[1:2] +sameVersion <- version12 == version_local12 +if (any(sameVersion == FALSE)) { +# generate message +version_msg <- strsplit(strsplit(page[versionLine], split = '

')[[1]][2], split = '

')[[1]] +infoLine_start <- versionLine + 1 +infoLine_end <- grep('

For historical releases and the introduction of updates about each version', page) +info_msg <- '' +for (infoLine in infoLine_start:infoLine_end) { +info_msg <- c(info_msg, strsplit(strsplit(page[infoLine], split = '

')[[1]][2], split = '

')[[1]]) +} +install_msg <- 'You can update by type in: devtools::install_gihub("Yuanchao-Xu/hyfo")' +message_out <- paste(version_msg, info_msg, install_msg, sep = '\n') +} else message_out <- NULL +return(message_out) +} +.onAttach <- function(libname, pkgname) { +message_out <- suppressWarnings(try(hyfoUpdates(), silent = TRUE)) +if (!is.null(message_out)) { +if (grepl('Version', message_out)) { +packageStartupMessage(message_out) +} +} +} +debug(hyfoUpdates) +hyfoUpdates() debug(hyfoUpdates) hyfoUpdates() -version12 -version_local12 sameVersion +version_msg +infoLine_start +infoLine_end +infoLine +strsplit(strsplit(page[infoLine], split = '

')[[1]][2], split = '

')[[1]] +strsplit(strsplit(page[57], split = '

')[[1]][2], split = '

')[[1]] +strsplit(strsplit(page[58], split = '

')[[1]][2], split = '

')[[1]] +strsplit(strsplit(page[59], split = '

')[[1]][2], split = '

')[[1]] +page[59] +strsplit(page[59], split = '

')[[1]][2], split = '

' +strsplit(page[59], split = '

')[[1]][2], split = '

') +strsplit(page[59], split = '

')[[1]][2] +strsplit(page[59], split = '

') +page[59] +page[58] +page[57] +page[56:64] ## For package updates information #' @importFrom utils packageDescription hyfoUpdates <- function(){ @@ -482,15 +240,95 @@ version_local <- packageDescription("hyfo")$Version version12 <- unlist(strsplit(version, split = "[.]"))[1:2] version_local12 <- unlist(strsplit(version_local, split = "[.]"))[1:2] sameVersion <- version12 == version_local12 +if (any(sameVersion == FALSE)) { # generate message version_msg <- strsplit(strsplit(page[versionLine], split = '

')[[1]][2], split = '

')[[1]] -infoLine <- versionLine + 2 -info_msg <- strsplit(strsplit(page[infoLine], split = '

')[[1]][2], split = '

')[[1]] +infoLine_start <- versionLine + 1 +infoLine_end <- grep('

For historical releases and the introduction of updates about each version', page) +info_msg <- '' +for (infoLine in infoLine_start:infoLine_end) { +info_msg <- c(info_msg, strsplit(strsplit(page[infoLine], split = '<*>')[[1]][2], split = '

')[[1]]) +} install_msg <- 'You can update by type in: devtools::install_gihub("Yuanchao-Xu/hyfo")' -message_out <- NULL +message_out <- paste(version_msg, info_msg, install_msg, sep = '\n') +} else message_out <- NULL +return(message_out) +} +.onAttach <- function(libname, pkgname) { +message_out <- suppressWarnings(try(hyfoUpdates(), silent = TRUE)) +if (!is.null(message_out)) { +if (grepl('Version', message_out)) { +packageStartupMessage(message_out) +} +} +} +## For package updates information +#' @importFrom utils packageDescription +hyfoUpdates <- function(){ +page <- readLines('http://yuanchao-xu.github.io/hyfo/') +updatesLine <- grep('id=\\"updates"', page) +versionLine <- updatesLine + 2 +version <- unlist(strsplit(page[versionLine], split = ' '))[2] +version_local <- packageDescription("hyfo")$Version +# the first tow digit is the most important part of the version +version12 <- unlist(strsplit(version, split = "[.]"))[1:2] +version_local12 <- unlist(strsplit(version_local, split = "[.]"))[1:2] +sameVersion <- version12 == version_local12 if (any(sameVersion == FALSE)) { +# generate message +version_msg <- strsplit(strsplit(page[versionLine], split = '

')[[1]][2], split = '

')[[1]] +infoLine_start <- versionLine + 1 +infoLine_end <- grep('

For historical releases and the introduction of updates about each version', page) +info_msg <- '' +for (infoLine in infoLine_start:infoLine_end) { +info_msg <- c(info_msg, strsplit(strsplit(page[infoLine], split = '<*>')[[1]][2], split = '

')[[1]]) +} +install_msg <- 'You can update by type in: devtools::install_gihub("Yuanchao-Xu/hyfo")' message_out <- paste(version_msg, info_msg, install_msg, sep = '\n') +} else message_out <- NULL +return(message_out) +} +.onAttach <- function(libname, pkgname) { +message_out <- suppressWarnings(try(hyfoUpdates(), silent = TRUE)) +if (!is.null(message_out)) { +if (grepl('Version', message_out)) { +packageStartupMessage(message_out) } +} +} +debug(hyfoUpdates) +hyfoUpdates() +strsplit(strsplit(page[infoLine], split = '<*>')[[1]][2], split = '<*>')[[1]] +strsplit(strsplit(page[infoLine], split = '<*>')[[1]][2] +strsplit(page[infoLine], split = '<*>')[[1]][2] +strsplit(page[infoLine], split = '>')[[1]][2] +strsplit(page[infoLine], split = '>') +strsplit(page[58], split = '>') +strsplit(page[infoLine], split = '>')Q +## For package updates information +#' @importFrom utils packageDescription +hyfoUpdates <- function(){ +page <- readLines('http://yuanchao-xu.github.io/hyfo/') +updatesLine <- grep('id=\\"updates"', page) +versionLine <- updatesLine + 2 +version <- unlist(strsplit(page[versionLine], split = ' '))[2] +version_local <- packageDescription("hyfo")$Version +# the first tow digit is the most important part of the version +version12 <- unlist(strsplit(version, split = "[.]"))[1:2] +version_local12 <- unlist(strsplit(version_local, split = "[.]"))[1:2] +sameVersion <- version12 == version_local12 +if (any(sameVersion == FALSE)) { +# generate message +version_msg <- strsplit(strsplit(page[versionLine], split = '

')[[1]][2], split = '

')[[1]] +infoLine_start <- versionLine + 2 +infoLine_end <- grep('

For historical releases and the introduction of updates about each version', page) - 1 +info_msg <- '' +for (infoLine in infoLine_start:infoLine_end) { +info_msg <- c(info_msg, strsplit(strsplit(page[infoLine], split = '>')[[1]][2], split = '<')[[1]]) +} +install_msg <- 'You can update by type in: devtools::install_gihub("Yuanchao-Xu/hyfo")' +message_out <- paste(version_msg, info_msg, install_msg, sep = '\n') +} else message_out <- NULL return(message_out) } .onAttach <- function(libname, pkgname) { @@ -503,10 +341,172 @@ packageStartupMessage(message_out) } debug(hyfoUpdates) hyfoUpdates() -N -sameVersion -message_out +infor_msg +info_msg +info_msg +strsplit(strsplit(page[infoLine], split = '>')[[1]][2], split = '<')[[1]] +strsplit(strsplit(page[infoLine], split = '>')[[1]][2], split = '<')[[1]][1] +q +## For package updates information +#' @importFrom utils packageDescription +hyfoUpdates <- function(){ +page <- readLines('http://yuanchao-xu.github.io/hyfo/') +updatesLine <- grep('id=\\"updates"', page) +versionLine <- updatesLine + 2 +version <- unlist(strsplit(page[versionLine], split = ' '))[2] +version_local <- packageDescription("hyfo")$Version +# the first tow digit is the most important part of the version +version12 <- unlist(strsplit(version, split = "[.]"))[1:2] +version_local12 <- unlist(strsplit(version_local, split = "[.]"))[1:2] +sameVersion <- version12 == version_local12 +if (any(sameVersion == FALSE)) { +# generate message +version_msg <- strsplit(strsplit(page[versionLine], split = '

')[[1]][2], split = '

')[[1]] +infoLine_start <- versionLine + 2 +infoLine_end <- grep('

For historical releases and the introduction of updates about each version', page) - 1 +info_msg <- '' +for (infoLine in infoLine_start:infoLine_end) { +info_line <- strsplit(strsplit(page[infoLine], split = '>')[[1]][2], split = '<')[[1]][1] +if (!is.na(info_line)) info_msg <- c(info_msg, info_line) +} +install_msg <- 'You can update by type in: devtools::install_gihub("Yuanchao-Xu/hyfo")' +message_out <- paste(version_msg, info_msg, install_msg, sep = '\n') +} else message_out <- NULL +return(message_out) +} +.onAttach <- function(libname, pkgname) { +message_out <- suppressWarnings(try(hyfoUpdates(), silent = TRUE)) +if (!is.null(message_out)) { +if (grepl('Version', message_out)) { +packageStartupMessage(message_out) +} +} +} +debug(hyfoUpdates) +hyfoUpdates() +info_msg +## For package updates information +#' @importFrom utils packageDescription +hyfoUpdates <- function(){ +page <- readLines('http://yuanchao-xu.github.io/hyfo/') +updatesLine <- grep('id=\\"updates"', page) +versionLine <- updatesLine + 2 +version <- unlist(strsplit(page[versionLine], split = ' '))[2] +version_local <- packageDescription("hyfo")$Version +# the first tow digit is the most important part of the version +version12 <- unlist(strsplit(version, split = "[.]"))[1:2] +version_local12 <- unlist(strsplit(version_local, split = "[.]"))[1:2] +sameVersion <- version12 == version_local12 +if (any(sameVersion == FALSE)) { +# generate message +version_msg <- strsplit(strsplit(page[versionLine], split = '

')[[1]][2], split = '

')[[1]] +infoLine_start <- versionLine + 2 +infoLine_end <- grep('

For historical releases and the introduction of updates about each version', page) - 1 +info_msg <- character() +for (infoLine in infoLine_start:infoLine_end) { +info_line <- strsplit(strsplit(page[infoLine], split = '>')[[1]][2], split = '<')[[1]][1] +if (!is.na(info_line)) info_msg <- c(info_msg, info_line) +} +install_msg <- 'You can update by type in: devtools::install_gihub("Yuanchao-Xu/hyfo")' +message_out <- paste(version_msg, info_msg, install_msg, sep = '\n') +} else message_out <- NULL +return(message_out) +} +.onAttach <- function(libname, pkgname) { +message_out <- suppressWarnings(try(hyfoUpdates(), silent = TRUE)) +if (!is.null(message_out)) { +if (grepl('Version', message_out)) { +packageStartupMessage(message_out) +} +} +} +debug(hyfoUpdates) +## For package updates information +#' @importFrom utils packageDescription +hyfoUpdates <- function(){ +page <- readLines('http://yuanchao-xu.github.io/hyfo/') +updatesLine <- grep('id=\\"updates"', page) +versionLine <- updatesLine + 2 +version <- unlist(strsplit(page[versionLine], split = ' '))[2] +version_local <- packageDescription("hyfo")$Version +# the first tow digit is the most important part of the version +version12 <- unlist(strsplit(version, split = "[.]"))[1:2] +version_local12 <- unlist(strsplit(version_local, split = "[.]"))[1:2] +sameVersion <- version12 == version_local12 +if (any(sameVersion == FALSE)) { +# generate message +version_msg <- strsplit(strsplit(page[versionLine], split = '

')[[1]][2], split = '

')[[1]] +infoLine_start <- versionLine + 2 +infoLine_end <- grep('

For historical releases and the introduction of updates about each version', page) - 1 +info_msg <- character() +for (infoLine in infoLine_start:infoLine_end) { +info_line <- strsplit(strsplit(page[infoLine], split = '>')[[1]][2], split = '<')[[1]][1] +if (!is.na(info_line)) info_msg <- c(info_msg, info_line) +} +install_msg <- 'You can update by type in: devtools::install_gihub("Yuanchao-Xu/hyfo")' +message_out <- paste(version_msg, info_msg, install_msg, sep = '\n') +} else message_out <- NULL +return(message_out) +} +.onAttach <- function(libname, pkgname) { +message_out <- suppressWarnings(try(hyfoUpdates(), silent = TRUE)) +if (!is.null(message_out)) { +if (grepl('Version', message_out)) { +packageStartupMessage(message_out) +} +} +} +debug(hyfoUpdates) +hyfoUpdates() +info_msg +devtools::document() +hyfoUpdates() +info_msg +version_msg +?paste +paste("1st", "2nd", "3rd", collapse = ", ") +paste("1st", "2nd", "3rd", sep = ", ") +paste("1st", "2nd", "3rd", collapse = ", dafa") +paste("1st", "2nd", "3rd", collapse = "dafa") +paste("1st", "2nd", "3rd", collapse = "") +paste("1st", "2nd", "3rd") +paste("1st", "2nd", "3rd", collapse = "-") +paste("1dafst", "2dasfnd", "dafa3rd", collapse = "-") +## For package updates information +#' @importFrom utils packageDescription +hyfoUpdates <- function(){ +page <- readLines('http://yuanchao-xu.github.io/hyfo/') +updatesLine <- grep('id=\\"updates"', page) +versionLine <- updatesLine + 2 +version <- unlist(strsplit(page[versionLine], split = ' '))[2] +version_local <- packageDescription("hyfo")$Version +# the first tow digit is the most important part of the version +version12 <- unlist(strsplit(version, split = "[.]"))[1:2] +version_local12 <- unlist(strsplit(version_local, split = "[.]"))[1:2] +sameVersion <- version12 == version_local12 +if (any(sameVersion == FALSE)) { +# generate message +version_msg <- strsplit(strsplit(page[versionLine], split = '

')[[1]][2], split = '

')[[1]] +infoLine_start <- versionLine + 2 +infoLine_end <- grep('

For historical releases and the introduction of updates about each version', page) - 1 +info_msg <- character() +for (infoLine in infoLine_start:infoLine_end) { +info_line <- strsplit(strsplit(page[infoLine], split = '>')[[1]][2], split = '<')[[1]][1] +if (!is.na(info_line)) info_msg <- c(info_msg, info_line) +} +install_msg <- 'You can update by type in: devtools::install_gihub("Yuanchao-Xu/hyfo")' +message_out <- paste(version_msg, paste(info_msg, collapse = '\n'), install_msg, sep = '\n') +} else message_out <- NULL +return(message_out) +} +.onAttach <- function(libname, pkgname) { +message_out <- suppressWarnings(try(hyfoUpdates(), silent = TRUE)) +if (!is.null(message_out)) { +if (grepl('Version', message_out)) { +packageStartupMessage(message_out) +} +} +} devtools::document() devtools::check() -devtools::check() -library(hyfo) +devtools::document() diff --git a/.Rproj.user/D53FD3E6/pcs/find-in-files.pper b/.Rproj.user/D53FD3E6/pcs/find-in-files.pper index 53dddf6..d1dda22 100644 --- a/.Rproj.user/D53FD3E6/pcs/find-in-files.pper +++ b/.Rproj.user/D53FD3E6/pcs/find-in-files.pper @@ -4,7 +4,7 @@ "filePatterns" : [ ], "path" : "E:/1/R/hyfo/R", - "query" : "default method is", + "query" : "paste", "regex" : false } } \ No newline at end of file diff --git a/.Rproj.user/D53FD3E6/pcs/source-pane.pper b/.Rproj.user/D53FD3E6/pcs/source-pane.pper index 3249574..70829f6 100644 --- a/.Rproj.user/D53FD3E6/pcs/source-pane.pper +++ b/.Rproj.user/D53FD3E6/pcs/source-pane.pper @@ -1,3 +1,3 @@ { - "activeTab" : 2 + "activeTab" : 1 } \ No newline at end of file diff --git a/.Rproj.user/D53FD3E6/pcs/workbench-pane.pper b/.Rproj.user/D53FD3E6/pcs/workbench-pane.pper index 7135db5..6785be6 100644 --- a/.Rproj.user/D53FD3E6/pcs/workbench-pane.pper +++ b/.Rproj.user/D53FD3E6/pcs/workbench-pane.pper @@ -1,4 +1,4 @@ { "TabSet1" : 0, - "TabSet2" : 0 + "TabSet2" : 3 } \ No newline at end of file diff --git a/.Rproj.user/D53FD3E6/sdb/per/t/222F1822 b/.Rproj.user/D53FD3E6/sdb/per/t/222F1822 index 4f5c521..7d95f90 100644 --- a/.Rproj.user/D53FD3E6/sdb/per/t/222F1822 +++ b/.Rproj.user/D53FD3E6/sdb/per/t/222F1822 @@ -1,12 +1,12 @@ { - "contents" : "hyfo 1.3.1\n==========\nDate: 2015.11.3\n\n- new generic function biasCorrect, extractPeriod, resample, getAnnual, getPreciBar added,\n No need to designate inputtype any more, R will detect automatically.\n- coordinates conversion function extracted.\n- new user manual about real time bias correction and resample added.\n\n\n\nhyfo 1.2.9\n==========\nDate: 2015.10.30\n\n- new biasFactor S4 class added, to avoid set the input type every time.\n- operational bias correction has been changed to generic function.\n- news file added.\n\n\n\nhyfo 1.2.8\n==========\nDate: 2015.10.10\n\n- operational bias correction added, in normal function.", + "contents" : "hyfo 1.3.1\n==========\nDate: 2015.11.3\n\n- new generic function biasCorrect, extractPeriod, resample, getAnnual, getPreciBar added. No need to designate inputtype any more, R will detect automatically.\n- coordinates conversion function extracted.\n- new user manual about real time bias correction and resample added.\n\n\n\nhyfo 1.2.9\n==========\nDate: 2015.10.30\n\n- new biasFactor S4 class added, to avoid set the input type every time.\n- operational bias correction has been changed to generic function.\n- news file added.\n\n\n\nhyfo 1.2.8\n==========\nDate: 2015.10.10\n\n- operational bias correction added, in normal function.", "created" : 1446423165783.000, "dirty" : false, "encoding" : "ASCII", "folds" : "", - "hash" : "884000713", + "hash" : "2337745590", "id" : "222F1822", - "lastKnownWriteTime" : 1446513924, + "lastKnownWriteTime" : 1446560725, "path" : "E:/1/R/hyfo/NEWS", "project_path" : "NEWS", "properties" : { diff --git a/.Rproj.user/D53FD3E6/sdb/per/t/3C7CEB0F b/.Rproj.user/D53FD3E6/sdb/per/t/3C7CEB0F deleted file mode 100644 index 91a4179..0000000 --- a/.Rproj.user/D53FD3E6/sdb/per/t/3C7CEB0F +++ /dev/null @@ -1,17 +0,0 @@ -{ - "contents" : "\n\n\n\n#' Get bias factor for multi/operational/real time bias correction.\n#' \n#' When you do multi/operational/real time bias correction. It's too expensive\n#' to input hindcast and obs every time. Especially when you have a long period of hindcast\n#' and obs, but only a short period of frc, it's too unecessary to read and compute hindcast\n#' and obs everytime. Therefore, biasFactor is designed. Using \\code{getBiasFactor}, you can\n#' get the biasFactor with hindcast and observation, then you can use \\code{applyBiasFactor} to \n#' apply the biasFactor to different forecasts. \n#' \n#' \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 method 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, default is 'multi'. More info check \n#' details.\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#' \n#' @seealso \\code{\\link{biasCorrect}} for method used in bias correction.\n#' \\code{\\link{applyBiasFactor}}, for the second part.\n#' \n#' @details \n#' \n#' Information about the method and how biasCorrect works can be found in \\code{\\link{biasCorrect}}\n#' \n#' \\strong{why use biasFactor}\n#' \n#' As for forecasting, for daily data, there is usually no need to have\n#' different bias factor every different day. You can calculate one bisa factor using a long\n#' period of hindcast and obs, and apply that factor to different frc.\n#' \n#' For example,\n#' \n#' You have 10 years of hindcast and observation. you want to do bias correction for some \n#' forecasting product, e.g. system 4. For system 4, each month, you will get a new forecast\n#' about the future 6 months. So if you want to do the real time bias correction, you have to\n#' take the 10 years of hindcast and observation data with you, and run \\code{biasCorrect} every\n#' time you get a new forecast. That's too expensive.\n#' \n#' For some practical use in forecasting, there isn't a so high demand for accuracy. E.g.,\n#' Maybe for February and March, you can use the same biasFactor, no need to do the computation \n#' again. \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#' \n#' # Then we will use nc data as forecasting data, and use itself as hindcast data,\n#' # use tgridData as observation.\n#' \n#' biasFactor <- getBiasFactor(nc, tgridData)\n#' newFrc <- applyBiasFactor(nc, biasFactor)\n#' \n#' biasFactor <- getBiasFactor(nc, tgridData, method = 'eqm', extrapolate = 'constant',\n#' preci = TRUE)\n#' # This method needs obs input.\n#' newFrc <- applyBiasFactor(nc, biasFactor, obs = tgridData)\n#' \n#' biasFactor <- getBiasFactor(nc, tgridData, method = 'gqm', preci = TRUE)\n#' newFrc <- applyBiasFactor(nc, biasFactor) \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\n#' biasFactor <- getBiasFactor(hindcast, obs)\n#' frc_new <- applyBiasFactor(frc, biasFactor)\n#' \n#' # for precipitation data, extra process needs to be executed, so you have to tell\n#' # the program to it is a precipitation data.\n#' \n#' biasFactor <- getBiasFactor(hindcast, obs, preci = TRUE)\n#' frc_new1 <- applyBiasFactor(frc, biasFactor)\n#' \n#' # You can use other methods to biascorrect, e.g. delta method. \n#' biasFactor <- getBiasFactor(hindcast, obs, method = 'delta')\n#' # delta method needs obs input.\n#' frc_new2 <- applyBiasFactor(frc, biasFactor, obs = obs)\n#' \n#' # \n#' biasFactor <- getBiasFactor(hindcast, obs, method = 'eqm', preci = TRUE)\n#' # eqm needs obs input\n#' frc_new3 <- applyBiasFactor(frc, biasFactor, obs = obs)\n#' \n#' biasFactor <- getBiasFactor(hindcast, obs, method = 'gqm', preci = TRUE)\n#' frc_new4 <- applyBiasFactor(frc, biasFactor)\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://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 }, S. Herrera \\email{sixto@@predictia.es }\n#' \n#' @importFrom methods setMethod\n#' @export\n#' \n#' \n# debug by trace(\"getBiasFactor\", browser, exit=browser, signature = c(\"list\", \"list\"))\nsetGeneric('getBiasFactor', function(hindcast, obs, method = 'scaling', scaleType = 'multi', \n preci = FALSE, prThreshold = 0, extrapolate = 'no') {\n standardGeneric('getBiasFactor')\n})\n\n#' @describeIn getBiasFactor\nsetMethod('getBiasFactor', signature('data.frame', 'data.frame'), \n function(hindcast, obs, method, scaleType, preci, prThreshold, extrapolate) {\n result <- getBiasFactor.TS(hindcast, obs, method, scaleType, preci, prThreshold, extrapolate)\n return(result)\n })\n\n\n# This is for the grid file from downscaleR\n#' @describeIn getBiasFactor\n#' @importFrom methods new\nsetMethod('getBiasFactor', signature('list', 'list'), \n function(hindcast, obs, method, scaleType, preci, prThreshold, extrapolate) {\n result <- getBiasFactor.list(hindcast, obs, method, scaleType, preci, prThreshold, extrapolate)\n return(result)\n })\n\n\n\n\n#' Apply bias factor to different forecasts for multi/operational/real time bias correction.\n#' \n#' When you do multi/operational/real time bias correction. It's too expensive\n#' to input hindcast and obs every time. Especially when you have a long period of hindcast\n#' and obs, but only a short period of frc, it's too unecessary to read and compute hindcast\n#' and obs everytime. Therefore, biasFactor is designed. Using \\code{getBiasFactor}, you can\n#' get the biasFactor with hindcast and observation, then you can use \\code{applyBiasFactor} to \n#' apply the biasFactor to different forecasts. \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 frc data. Check details for more information.\n#' @param biasFactor a file containing all the information of the calibration, will be\n#' applied to different forecasts.\n#' @param obs for some methods, observation input is necessary. obs is 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. Default value is NULL.\n#' @seealso \\code{\\link{biasCorrect}} for method used in bias correction. \n#' \\code{\\link{getBiasFactor}}, for the first part.\n#' \n#' @details \n#' \n#' Information about the method and how biasCorrect works can be found in \\code{\\link{biasCorrect}}\n#' \n#' \\strong{why use biasFactor}\n#' \n#' As for forecasting, for daily data, there is usually no need to have\n#' different bias factor every different day. You can calculate one bisa factor using a long\n#' period of hindcast and obs, and apply that factor to different frc.\n#' \n#' For example,\n#' \n#' You have 10 years of hindcast and observation. you want to do bias correction for some \n#' forecasting product, e.g. system 4. For system 4, each month, you will get a new forecast\n#' about the future 6 months. So if you want to do the real time bias correction, you have to\n#' take the 10 years of hindcast and observation data with you, and run \\code{biasCorrect} every\n#' time you get a new forecast. That's too expensive.\n#' \n#' For some practical use in forecasting, there isn't a so high demand for accuracy. E.g.,\n#' Maybe for February and March, you can use the same biasFactor, no need to do the computation \n#' again. \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#' \n#' # Then we will use nc data as forecasting data, and use itself as hindcast data,\n#' # use tgridData as observation.\n#' \n#' biasFactor <- getBiasFactor(nc, tgridData)\n#' newFrc <- applyBiasFactor(nc, biasFactor)\n#' \n#' biasFactor <- getBiasFactor(nc, tgridData, method = 'eqm', extrapolate = 'constant',\n#' preci = TRUE)\n#' # This method needs obs input.\n#' newFrc <- applyBiasFactor(nc, biasFactor, obs = tgridData)\n#' \n#' biasFactor <- getBiasFactor(nc, tgridData, method = 'gqm', preci = TRUE)\n#' newFrc <- applyBiasFactor(nc, biasFactor) \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\n#' biasFactor <- getBiasFactor(hindcast, obs)\n#' frc_new <- applyBiasFactor(frc, biasFactor)\n#' \n#' # for precipitation data, extra process needs to be executed, so you have to tell\n#' # the program to it is a precipitation data.\n#' \n#' biasFactor <- getBiasFactor(hindcast, obs, preci = TRUE)\n#' frc_new1 <- applyBiasFactor(frc, biasFactor)\n#' \n#' # You can use other methods to biascorrect, e.g. delta method. \n#' biasFactor <- getBiasFactor(hindcast, obs, method = 'delta')\n#' # delta method needs obs input.\n#' frc_new2 <- applyBiasFactor(frc, biasFactor, obs = obs)\n#' \n#' # \n#' biasFactor <- getBiasFactor(hindcast, obs, method = 'eqm', preci = TRUE)\n#' # eqm needs obs input\n#' frc_new3 <- applyBiasFactor(frc, biasFactor, obs = obs)\n#' \n#' biasFactor <- getBiasFactor(hindcast, obs, method = 'gqm', preci = TRUE)\n#' frc_new4 <- applyBiasFactor(frc, biasFactor)\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://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 }, S. Herrera \\email{sixto@@predictia.es }\n#' \n#' @export\nsetGeneric('applyBiasFactor', function(frc, biasFactor, obs = NULL) {\n standardGeneric('applyBiasFactor')\n})\n\n#' @describeIn applyBiasFactor\n#' @importFrom methods setMethod\nsetMethod('applyBiasFactor', signature('data.frame', 'biasFactor'), \n function(frc, biasFactor, obs) {\n result <- applyBiasFactor.TS(frc, biasFactor, obs)\n return(result)\n })\n \n#' @describeIn applyBiasFactor\n#' @importFrom methods setMethod\nsetMethod('applyBiasFactor', signature('list', 'biasFactor.hyfo'), \n function(frc, biasFactor, obs) {\n result <- applyBiasFactor.list(frc, biasFactor, obs)\n return(result)\n })\n\n\n### generic functions\ngetBiasFactor.TS <- function(hindcast, obs, method, scaleType, preci, prThreshold, extrapolate) {\n \n if (!grepl('-|/', obs[1, 1]) | !grepl('-|/', hindcast[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 \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(hindcast)\n \n # For every column, it's biascorrected respectively.\n biasFactor <- lapply(2:n, function(x) getBiasFactor_core(hindcast[, x], obs[, 2], method = method,\n scaleType = scaleType, preci = preci, prThreshold = prThreshold, \n extrapolate = extrapolate))\n if (n - 1 > 1) {\n biasFactor_all <- new('biasFactor.multiMember', biasFactor = biasFactor, memberDim = n - 1,\n method = method, preci = preci, prThreshold = prThreshold, scaleType = scaleType, \n extrapolate = extrapolate)\n \n } else {\n biasFactor_all <- new('biasFactor', biasFactor = biasFactor, method = method, \n preci = preci, prThreshold = prThreshold, scaleType = scaleType, \n extrapolate = extrapolate)\n }\n \n return(biasFactor_all)\n}\n\ngetBiasFactor.list <- function(hindcast, obs, method, scaleType, preci, prThreshold, extrapolate) {\n \n ## Check if the data is a hyfo grid data.\n checkHyfo(hindcast, obs)\n \n hindcastData <- hindcast$Data\n obsData <- obs$Data\n \n ## save frc dimension order, at last, set the dimension back to original dimension\n hindcastDim <- attributes(hindcastData)$dimensions\n \n ## ajust the dimension into general dimension order.\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(hindcastData, obsData, dim = c('lon', 'lat'))\n \n \n # Now real bias correction is executed.\n \n memberIndex <- match('member', attributes(hindcastData)$dimensions)\n \n # For dataset that has a member part \n if (!is.na(memberIndex)) {\n \n hindcastData <- adjustDim(hindcastData, 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 biasFactor_all <- vector(mode = \"list\", length = dim(hindcastData)[4])\n for (member in 1:dim(hindcastData)[4]) {\n biasFactor_all[[member]] <- vector(mode = 'list', length = dim(hindcastData)[1])\n for (lon in 1:dim(hindcastData)[1]) {\n biasFactor_all[[member]][[lon]] <- vector(mode = 'list', length = dim(hindcastData)[2])\n for (lat in 1:dim(hindcastData)[2]) {\n biasFactor_all[[member]][[lon]][[lat]] <- getBiasFactor_core(hindcastData[lon, lat,, member], obsData[lon, lat,], method = method,\n scaleType = scaleType, preci = preci, prThreshold = prThreshold, \n extrapolate = extrapolate)\n }\n }\n }\n \n biasFactor <- new('biasFactor.hyfo', biasFactor = biasFactor_all, method = method, preci = preci,\n prThreshold = prThreshold, scaleType = scaleType, extrapolate = extrapolate, \n lonLatDim = calcuDim(hindcastData, dim = c('lon', 'lat')),\n memberDim = calcuDim(hindcastData, dim = 'member'))\n } else {\n \n hindcastData <- adjustDim(hindcastData, ref = c('lon', 'lat', 'time'))\n \n biasFactor_all <- vector(mode = 'list', length = dim(hindcastData)[1])\n for (lon in 1:dim(hindcastData)[1]) {\n biasFactor_all[[lon]] <- vector(mode = 'list', length = dim(hindcastData)[2]) \n for (lat in 1:dim(hindcastData)[2]) {\n biasFactor_all[[lon]][[lat]] <- getBiasFactor_core(hindcastData[lon, lat,], obsData[lon, lat,], method = method,\n scaleType = scaleType, preci = preci, prThreshold = prThreshold, \n extrapolate = extrapolate)\n }\n }\n biasFactor <- new('biasFactor.hyfo', biasFactor = biasFactor_all, method = method, preci = preci,\n prThreshold = prThreshold, scaleType = scaleType, extrapolate = extrapolate, \n lonLatDim = calcuDim(hindcastData, dim = c('lon', 'lat')))\n \n }\n \n return(biasFactor)\n}\n\napplyBiasFactor.TS <- function(frc, biasFactor, obs) {\n method <- biasFactor@method\n preci <- biasFactor@preci\n prThreshold <- biasFactor@prThreshold\n scaleType <- biasFactor@scaleType\n extrapolate <- biasFactor@extrapolate\n memberDim <- biasFactor@memberDim\n biasFactor <- biasFactor@biasFactor\n \n \n # First check if the first column is Date\n if (!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 # In this case more than one value columns exist in the dataset, both frc and hindcast.\n \n n <- ncol(frc)\n if (n-1 != memberDim) stop('frc and biasFactor have different members.')\n \n \n # For every column, it's biascorrected respectively.\n frc_data <- lapply(2:n, function(x) applyBiasFactor_core(frc[, x], biasFactor = biasFactor[[x - 1]], method = method,\n scaleType = scaleType, preci = preci, prThreshold = prThreshold, \n extrapolate = extrapolate, obs = obs[, 2]))\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 \n return(frc_new)\n \n}\n\napplyBiasFactor.list <- function(frc, biasFactor, obs) {\n method <- biasFactor@method\n preci <- biasFactor@preci\n prThreshold <- biasFactor@prThreshold\n scaleType <- biasFactor@scaleType\n extrapolate <- biasFactor@extrapolate\n lonLatDim <- biasFactor@lonLatDim\n memberDim <- biasFactor@memberDim\n biasFactor <- biasFactor@biasFactor\n \n ## Check if the data is a hyfo grid data.\n checkHyfo(frc)\n \n \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 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 if (!identical(calcuDim(frcData, dim = c('lon', 'lat')), lonLatDim)) {\n stop('frc data has different lon and lat from hindcast data.')\n }\n \n \n # Now real bias correction is executed.\n \n memberIndex <- match('member', attributes(frcData)$dimensions)\n \n # For dataset that has a member part \n if (!is.na(memberIndex)) {\n # check if frcData and hindcastData has the same dimension and length.\n if (calcuDim(frcData, dim = 'member') != memberDim) {\n stop('frc data has different member number from hindcast.')\n } \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] <- applyBiasFactor_core(frcData[lon, lat,,member], biasFactor = biasFactor[[member]][[lon]][[lat]], method = method,\n scaleType = scaleType, preci = preci, prThreshold = prThreshold, \n extrapolate = extrapolate, obs = obsData[lon, lat,])\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,] <- applyBiasFactor_core(frcData[lon, lat,], biasFactor = biasFactor[[lon]][[lat]], method = method,\n scaleType = scaleType, preci = preci, prThreshold = prThreshold, \n extrapolate = extrapolate, obs = obsData[lon, lat,])\n }\n }\n }\n \n frcData <- adjustDim(frcData, ref = frcDim)\n frc$Data <- frcData\n frc$biasCorrected_by <- method\n frc_new <- frc\n \n return(frc_new)\n}\n\n\n#################\n################# core functions for multi bias correction.\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#' }\ngetBiasFactor_core <- function(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(hindcast))) {\n warning('In this cell, hindcast or obs data is missing. No biasCorrection for this cell.')\n return(NA)\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 biasFactor <- getBiasFactor_core_delta(hindcast)\n } else if (method == 'scaling') {\n biasFactor <- getBiasFactor_core_scaling(hindcast, obs, scaleType)\n } else if (method == 'eqm') {\n # In this method, the value is bounded by the observation\n # Preci or not both have the same biasFactor\n if (preci == FALSE) {\n biasFactor <- getBiasFactor_core_eqm_nonPreci(hindcast, obs, extrapolate)\n } else {\n biasFactor <- getBiasFactor_core_eqm_preci(hindcast, obs, minHindcastPreci, extrapolate, prThreshold)\n }\n \n \n } else if (method == 'gqm') {\n if (preci == FALSE) stop ('gqm method only applys to precipitation, please set preci = T')\n biasFactor <- getBiasFactor_core_gqm(hindcast, obs, prThreshold, minHindcastPreci)\n }\n \n if (preci == TRUE) biasFactor$minHindcastPreci <- minHindcastPreci\n \n return(biasFactor)\n}\n\n\napplyBiasFactor_core <- function(frc, biasFactor, method, preci, prThreshold, scaleType,\n extrapolate, obs = NULL) {\n \n if (!any(!is.na(biasFactor))) {\n warning('In this cell, biasFactor is missing.No biasCorrection for this cell.')\n # here return NA or return the unprocessed frc, both are OK. But return NA is more\n # obvious for user.\n return(NA)\n }\n \n if (method == 'delta') {\n if (is.null(obs)) stop('This method needs obs input.')\n if (length(frc) != length(obs)) stop('This method needs frc data have the same length as obs data.')\n frc <- applyBiasFactor_core_delta(frc = frc, biasFactor = biasFactor, obs = obs)\n } else if (method == 'scaling') {\n frc <- applyBiasFactor_core_scaling(frc = frc, biasFactor = biasFactor, scaleType = scaleType)\n } else if (method == 'eqm') {\n if (is.null(obs)) stop('This method needs obs input.')\n if (preci == FALSE) {\n frc <- applyBiasFactor_core_eqm_nonPreci(frc = frc, biasFactor = biasFactor, extrapolate = extrapolate, \n obs = obs)\n } else {\n frc <- applyBiasFactor_core_eqm_preci(frc = frc, biasFactor = biasFactor, extrapolate = extrapolate, \n prThreshold = prThreshold, obs = obs)\n }\n } else if (method == 'gqm') {\n frc <- applyBiasFactor_core_gqm(frc = frc, biasFactor = biasFactor)\n }\n \n return(frc)\n}\n\n\ngetBiasFactor_core_delta <- function(hindcast) {\n biasFactor <- list()\n biasFactor$hindcastMean <- mean(hindcast, na.rm = TRUE)\n return(biasFactor)\n}\napplyBiasFactor_core_delta <- function(frc, biasFactor, obs) {\n hindcastMean <- biasFactor$hindcastMean\n frcMean <- mean(frc, na.rm = TRUE)\n return(obs - hindcastMean + frcMean)\n}\n\ngetBiasFactor_core_scaling <- function(hindcast, obs, scaleType) {\n biasFactor <- list()\n \n hindcastMean <- mean(hindcast, na.rm = TRUE)\n obsMean <- mean(obs, na.rm = TRUE)\n \n if (scaleType == 'multi') {\n biasFactor$scale <- obsMean / hindcastMean\n \n } else if (scaleType == 'add') {\n biasFactor$scale <- obsMean - hindcastMean\n }\n \n return(biasFactor)\n}\n\napplyBiasFactor_core_scaling <- function(frc, biasFactor, scaleType) {\n \n if (scaleType == 'multi') {\n frc <- frc * biasFactor$scale\n \n } else if (scaleType == 'add') {\n frc <- frc + biasFactor$scale\n }\n return(frc)\n}\n\ngetBiasFactor_core_eqm_nonPreci <- function(hindcast, obs, extrapolate) {\n \n biasFactor <- list()\n biasFactor$ecdfHindcast <- ecdf(hindcast)\n \n if (extrapolate == 'constant') {\n biasFactor$maxHindcast <- max(hindcast, na.rm = TRUE)\n biasFactor$minHindcast <- min(hindcast, na.rm = TRUE)\n biasFactor$higherIndex_dif <- biasFactor$maxHindcast - max(obs, na.rm = TRUE)\n biasFactor$lowerIndex_dif <- biasFactor$minHindcast - min(obs, na.rm = TRUE)\n }\n \n return(biasFactor)\n}\n\ngetBiasFactor_core_eqm_preci <- function(hindcast, obs, minHindcastPreci, extrapolate,\n prThreshold) {\n \n biasFactor <- list()\n biasFactor$ecdfHindcast <- ecdf(hindcast[hindcast > minHindcastPreci])\n \n if (extrapolate == 'constant') {\n biasFactor$maxHindcast <- max(hindcast, na.rm = TRUE)\n biasFactor$minHindcast <- min(hindcast, na.rm = TRUE)\n biasFactor$higherIndex_dif <- biasFactor$maxHindcast - max(obs, na.rm = TRUE)\n biasFactor$lowerIndex_dif <- biasFactor$minHindcast - min(obs, nna.rm = TRUE)\n }\n biasFactor$availableHindcastLength <- length(which(hindcast > minHindcastPreci))\n \n # drizzle parameter 1\n biasFactor$drizzleP1 <- min(hindcast[hindcast > minHindcastPreci], na.rm = TRUE)\n # biasFactor$prThreshold <- prThreshold\n return(biasFactor)\n}\n\napplyBiasFactor_core_eqm_nonPreci <- function(frc, biasFactor, extrapolate, obs) {\n ecdfHindcast <- biasFactor$ecdfHindcast\n \n if (extrapolate == 'constant') {\n higherIndex <- which(frc > biasFactor$maxHindcast)\n lowerIndex <- which(frc < biasFactor$minHindcast)\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 \n frc[higherIndex] <- frc[higherIndex] - biasFactor$higherIndex_dif\n }\n \n if (length(lowerIndex) > 0) {\n \n frc[lowerIndex] <- frc[lowerIndex] - biasFactor$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\n#' @importFrom stats quantile\napplyBiasFactor_core_eqm_preci <- function(frc, biasFactor, extrapolate, prThreshold, obs) {\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 (biasFactor$availableHindcastLength > 0) {\n \n ecdfHindcast <- biasFactor$ecdfHindcast\n \n noRain <- which(frc <= biasFactor$minHindcastPreci & !is.na(frc))\n rain <- which(frc > biasFactor$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 > biasFactor$minHindcastPreci & frc <= biasFactor$drizzleP1 & !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] > biasFactor$maxHindcast)\n lowerIndex <- which(frc[rain] < biasFactor$minHindcast)\n \n extrapolateIndex <- c(higherIndex, lowerIndex)\n non_extrapolateIndex <- setdiff(1:length(rain), extrapolateIndex)\n \n if (length(higherIndex) > 0) {\n frc[rain[higherIndex]] <- frc[higherIndex] - biasFactor$higherIndex_dif\n }\n \n if (length(lowerIndex) > 0) {\n frc[rain[lowerIndex]] <- frc[lowerIndex] - biasFactor$lowerIndex_dif\n }\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 \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 > biasFactor$drizzleP1 & !is.na(frc))], \n 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 <= biasFactor$minHindcastPreci & !is.na(frc))\n rain <- which(frc > biasFactor$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\n#' @importFrom MASS fitdistr\ngetBiasFactor_core_gqm <- function(hindcast, obs, prThreshold, minHindcastPreci) {\n if (any(obs > prThreshold)) {\n biasFactor <- list()\n ind <- which(obs > prThreshold & !is.na(obs))\n obsGamma <- fitdistr(obs[ind],\"gamma\")\n biasFactor$obsShape <- obsGamma$estimate[1]\n biasFactor$obsRate <- obsGamma$estimate[2]\n \n ind <- which(hindcast > 0 & !is.na(hindcast))\n hindcastGamma <- fitdistr(hindcast[ind],\"gamma\")\n biasFactor$hindcastShape <- hindcastGamma$estimate[1]\n biasFactor$hindcastRate <- hindcastGamma$estimate[2]\n biasFactor$minHindcastPreci <- minHindcastPreci\n \n } else {\n warning('All the observations of this cell(station) are lower than the threshold, \n no biasFactor returned.')\n biasFactor <- NA\n }\n return(biasFactor)\n}\n\n#' @importFrom stats pgamma qgamma\napplyBiasFactor_core_gqm <- function(frc, biasFactor) {\n \n rain <- which(frc > biasFactor$minHindcastPreci & !is.na(frc))\n noRain <- which(frc <= biasFactor$minHindcastPreci & !is.na(frc))\n \n probF <- pgamma(frc[rain], biasFactor$hindcastShape, rate = biasFactor$hindcastRate)\n frc[rain] <- qgamma(probF, biasFactor$obsShape, rate = biasFactor$obsRate)\n frc[noRain] <- 0\n \n return(frc)\n}", - "created" : 1446554862150.000, - "dirty" : false, - "encoding" : "ASCII", - "folds" : "", - "hash" : "2456348229", - "id" : "3C7CEB0F", - "lastKnownWriteTime" : 1446513731, - "path" : "E:/1/R/hyfo/R/multi-biasCorrect(generic).R", - "project_path" : "R/multi-biasCorrect(generic).R", - "properties" : { - }, - "relative_order" : 8, - "source_on_save" : false, - "type" : "r_source" -} \ No newline at end of file diff --git a/.Rproj.user/D53FD3E6/sdb/per/t/44AE41D2 b/.Rproj.user/D53FD3E6/sdb/per/t/44AE41D2 index 87d9c45..043f6aa 100644 --- a/.Rproj.user/D53FD3E6/sdb/per/t/44AE41D2 +++ b/.Rproj.user/D53FD3E6/sdb/per/t/44AE41D2 @@ -6,7 +6,7 @@ "folds" : "", "hash" : "3174204225", "id" : "44AE41D2", - "lastKnownWriteTime" : 1446556010, + "lastKnownWriteTime" : 1446562168, "path" : "E:/1/R/hyfo/NAMESPACE", "project_path" : "NAMESPACE", "properties" : { diff --git a/.Rproj.user/D53FD3E6/sdb/per/t/4B970497 b/.Rproj.user/D53FD3E6/sdb/per/t/4B970497 deleted file mode 100644 index 287205b..0000000 --- a/.Rproj.user/D53FD3E6/sdb/per/t/4B970497 +++ /dev/null @@ -1,17 +0,0 @@ -{ - "contents" : "#' Get mean rainfall data.\n#' \n#' Get mean rainfall data, e.g. mean annual rainfall, mean monthly rainfall and mean winter rainfall.\n#' \n#' @param inputTS A time series with only data column (1 column).\n#' @param method A string showing the method used to calculate mean value, e.g., \"annual\".\n#' more information please refer to details.\n#' @param yearIndex A NUMERIC ARRAY showing the year index of the time series.\n#' @param monthIndex A NUMERIC ARRAY showing the month index of the time series.\n#' @param fullResults A boolean showing whether the full results are shown, default is FALSE. If \n#' FALSE, only mean value will be returned, if TRUE, the sequence of values will be returned.\n#' @param omitNA A boolean showing in the calculation, whether NA is omitted, default is FALSE.\n#' @param plot A boolean showing whether the results will be plotted.\n#' @param ..., \\code{title, x, y} showing the title and x and y axis of the plot, shoud be a string.\n#' @details\n#' There are following methods to be selected, \n#' \"annual\": annual rainfall of each year is plotted. \n#' \"winter\", \"spring\", \"autumn\", \"summer\": seasonal rainfall of each year is plotted.\n#' Month(number 1 to 12): month rainfall of each year is plotted, e.g. march rainfall of each year.\n#' \"meanMonthly\": the mean monthly rainfall of each month over the whole period.\n#' \n#' Since \"winter\" is a crossing year, 12, 1, 2, 12 is in former year, and 1, 2 are in latter year.\n#' so winter belongs to the latter year.\n#' \n#' @return The mean value of the input time series or the full results before calculating mean.\n#' \n# data(testdl)\n# TS <- testdl[[1]]\n# year = as.numeric(format(TS[, 1], '%Y'))\n# month = as.numeric(format(TS[, 1], '%m'))\n# \n# # Get the mean spring precipitation.\n# a <- getMeanPreci(TS[, 2], method = 'spring', yearIndex = year, monthIndex = month)\n# a\n# \n# # Get the series of spring precipitation, set fullResults = TRUE.\n# a <- getMeanPreci(TS[, 2], method = 'spring', yearIndex = year, monthIndex = month,\n# fullResults = TRUE)\n# a\n# \n# # If missing value is excluded, set omitNA = TRUE.\n# a <- getMeanPreci(TS[, 2], method = 'winter', yearIndex = year, monthIndex = month,\n# omitNA = TRUE, fullResults = TRUE)\n# a\n# \n# # Get special month precipitation, e.g. march.\n# a <- getMeanPreci(TS[, 2], method = 3, yearIndex = year, monthIndex = month,\n# fullResults = TRUE)\n# a\n# \n# # We can also get annual precipitation.\n# a <- getMeanPreci(TS[, 2], method = 'annual', yearIndex = year, monthIndex = month,\n# fullResults = TRUE)\n\ngetMeanPreci <- function(inputTS, method = NULL, yearIndex = NULL, monthIndex = NULL,\n fullResults = FALSE, omitNA = TRUE, plot = FALSE, ...) {\n # First check if all the records are NA.\n if (any(!is.na(inputTS))) {\n #converting daily preci to the wanted preci.\n if (method == 'annual') {\n ###yearIndex <- startTime$year + 1900\n annualPreci <- tapply(inputTS, INDEX = yearIndex, FUN = sum, na.rm = omitNA)#ggplot is able not to show NA, so choose TRUE\n if (fullResults == TRUE) output <- annualPreci else output <- mean(annualPreci, na.rm = TRUE)\n \n } else if (method == 'meanMonthly') {\n \n monthlypreci <- tapply(inputTS, INDEX = list(yearIndex, monthIndex), FUN = sum, na.rm = omitNA)\n meanMonthlyPreci <- apply(monthlypreci, MARGIN = 2, FUN = mean, na.rm = TRUE)\n \n if (fullResults == TRUE) output <- meanMonthlyPreci else output <- mean(meanMonthlyPreci, na.rm = TRUE)\n \n }else if (method == 'winter') {\n# #winter is the most tricky part, because it starts from Dec to Feb next year, it's a year-crossing season,\n# #so we have to make some changes to the monthIndex\n# #e.g.data from 1950.1.1 - 2008.3.31 if we want to calculate the mean winter preci, to calculate winter month\n# #December, we have to move the yearIndex one month forwards or two months backwards, to make 12,1,2 in one year \n# ###yearIndex <- startTime$year + 1900\n# ###monthIndex <- startTime$mon + 1\n# \n# #we move the yearIndex one month backwards\n# yearIndex_new <- c(yearIndex[32:length(yearIndex)], rep(tail(yearIndex, 1), 31))\n# \n# winterIndex <- which(monthIndex == 12 | monthIndex == 1 | monthIndex == 2)\n# winterYear <- yearIndex_new[winterIndex]#this index is used for calculation\n# \n# #because we don't have 1949.Dec, so the first winter is not intact, so first two months are elemenated\n# \n# startIndex <- length(which(winterYear == yearIndex[1])) + 1\n# winterOfLastYear <- length(which(winterYear == tail(yearIndex, 1)))\n# if (winterOfLastYear > 91) {\n# endIndex <- length(winterYear) - 31 #in case the data set ends at Dec.31\n# \n# } else if (winterOfLastYear < 90) { # incase the data ends at Jan 31\n# endIndex <- length(winterYear) - length(which(winterYear == tail(yearIndex, 1)))\n# \n# } else {\n# endIndex <- length(winterYear)\n# }\n# \n# inputTS <- inputTS[winterIndex][startIndex:endIndex]#needs two process with inputPreci, first, extract\n# #the winter preci, second, delete first two month of 1950\n# \n# winterYear <- winterYear[startIndex:endIndex]#needs one process, delete two months \n# seasonalPreci <- tapply(inputTS,INDEX = winterYear, FUN = sum, na.rm = omitNA)\n \n # upper part is the older method saved as backup.\n \n matrix <- tapply(inputTS, INDEX = list(monthIndex, yearIndex), FUN = sum, na.rm = omitNA)\n col <- colnames(matrix)\n dec <- matrix['12',] # extract December.\n dec <- c(NA, dec[1:length(dec) - 1]) # rearrange December order to push it to next year.\n names(dec) <- col\n matrix <- rbind(dec, matrix[rownames(matrix) != '12', ]) \n seasonalPreci <- apply(matrix, MARGIN = 2, function(x) sum(x[c('dec', '1', '2')]))\n \n if (fullResults == TRUE) output <- seasonalPreci else output <- mean(seasonalPreci, na.rm = TRUE) \n \n } else if (method == 'spring') {\n \n# springIndex <- which(monthIndex == 3 | monthIndex == 4 | monthIndex == 5)\n# springYear <- yearIndex[springIndex]\n# inputTS <- inputTS[springIndex]\n# seasonalPreci <- tapply(inputTS, INDEX = springYear, FUN = sum, na.rm = omitNA)\n \n \n matrix <- tapply(inputTS, INDEX = list(monthIndex, yearIndex), FUN = sum, na.rm = omitNA)\n seasonalPreci <- apply(matrix, MARGIN = 2, function(x) sum(x[c('3', '4', '5')]))\n \n if (fullResults == TRUE) output <- seasonalPreci else output <- mean(seasonalPreci, na.rm = TRUE)\n \n } else if (method == 'summer') {\n \n matrix <- tapply(inputTS, INDEX = list(monthIndex, yearIndex), FUN = sum, na.rm = omitNA)\n seasonalPreci <- apply(matrix, MARGIN = 2, function(x) sum(x[c('6', '7', '8')]))\n \n if (fullResults == TRUE) output <- seasonalPreci else output <- mean(seasonalPreci, na.rm = TRUE)\n \n } else if (method == 'autumn') {\n \n matrix <- tapply(inputTS, INDEX = list(monthIndex, yearIndex), FUN = sum, na.rm = omitNA)\n seasonalPreci <- apply(matrix, MARGIN = 2, function(x) sum(x[c('9', '10', '11')]))\n\n if (fullResults == TRUE) output <- seasonalPreci else output <- mean(seasonalPreci, na.rm = TRUE)\n \n } else if (is.numeric(method)) {\n \n month <- method\n \n #check if month exist \n e <- match(month, unique(monthIndex))\n if (is.na(e)) {\n e1 <- paste(unique(monthIndex), collapse = ',')\n m <- paste('No input month exists in the dataset, choose among', e1)\n stop(m)\n }\n \n monthlyPreci <- tapply(inputTS, INDEX = list(yearIndex, monthIndex), \n FUN = sum, na.rm = omitNA)[, toString(month)]\n \n if (fullResults == TRUE) output <- monthlyPreci else output <- mean(monthlyPreci, na.rm = TRUE)\n }\n \n } else {\n output <- NA\n }\n\n if (plot == TRUE) {\n a <- data.frame(Date = names(output), value = output)\n \n theme_set(theme_bw())\n mainLayer <- with(a, {\n ggplot(a) +\n geom_bar(aes(x = Date, y = value), stat = 'identity', fill = 'cyan') +\n labs(empty = NULL, ...) +#in order to pass \"...\", arguments shouldn't be empty.\n theme(plot.title = element_text(size = rel(1.3), face = 'bold'),\n axis.title.x = element_text(size = rel(1.2)),\n axis.title.y = element_text(size = rel(1.2))) + \n theme(axis.text.x = element_text(angle = 90, hjust = 1))\n })\n \n print (mainLayer)\n }\n \n return(output)\n}\n", - "created" : 1446424240843.000, - "dirty" : false, - "encoding" : "ASCII", - "folds" : "", - "hash" : "1620316423", - "id" : "4B970497", - "lastKnownWriteTime" : 1446467132, - "path" : "E:/1/R/hyfo/R/getMeanPreci.R", - "project_path" : "R/getMeanPreci.R", - "properties" : { - }, - "relative_order" : 4, - "source_on_save" : false, - "type" : "r_source" -} \ No newline at end of file diff --git a/.Rproj.user/D53FD3E6/sdb/per/t/5ECF1EB5 b/.Rproj.user/D53FD3E6/sdb/per/t/5ECF1EB5 deleted file mode 100644 index 7552011..0000000 --- a/.Rproj.user/D53FD3E6/sdb/per/t/5ECF1EB5 +++ /dev/null @@ -1,18 +0,0 @@ -{ - "contents" : "#' Monthly data to daily data and the reverse conversion.\n#' \n#' @param data a hyfo grid data or a time series, with first column date, and second column value. The date column should\n#' follow the format in \\code{as.Date}, i.e. seperate with \"-\" or \"/\". Check details for more information.\n#' @param method A string showing whether you want to change a daily data to monthly data or monthly\n#' data to daily data.e.g. \"mon2day\" and \"day2mon\".\n#' @details \n#' Note, when you want to change daily data to monthly data, a new date column will be generated,\n#' usually the date column will be the middle date of each month, 15th, or 16th. However, if your \n#' time series doesn't start from the beginning of a month or ends to the end of a month, e.g. \n#' from 1999-3-14 to 2008-2-2, the first and last generated date could be wrong. Not only the date, but also the data, because you are \n#' not calculating based on a intact month. \n#' @return converted time series.\n#' @examples\n#' # Daily to monthly\n#' data(testdl)\n#' TS <- testdl[[2]] # Get daily data\n#' str(TS)\n#' TS_new <- resample(TS, method = 'day2mon')\n#' \n#' # Monthly to daily\n#' TS <- data.frame(Date = seq(as.Date('1999-9-15'), length = 30, by = '1 month'), \n#' runif(30, 3, 10))\n#' TS_new <- resample(TS, method = 'mon2day')\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#' nc_new <- resample(nc, 'day2mon')\n#' \n#' \n#' # More examples can be found in the user manual on http://yuanchao-xu.github.io/hyfo/\n#' \n#' @export\n#' @importFrom stats aggregate\n#' @references \n#' \n#' \\itemize{\n#' \\item R Core Team (2015). R: A language and environment for statistical computing. R Foundation for\n#' Statistical Computing, Vienna, Austria. URL http://www.R-project.org/.\n#' }\n#' \nsetGeneric('resample', function(data, method) {\n standardGeneric('resample')\n})\n\n\n#' @describeIn resample\nsetMethod('resample', signature('data.frame'),\n function(data, method) {\n result <- resample.TS(data, method)\n return(result)\n })\n\n#' @describeIn resample\nsetMethod('resample', signature('list'),\n function(data, method) {\n result <- resample.list(data, method)\n return(result)\n })\n\n\n\n#' @importFrom stats aggregate\nresample.TS <- function(TS, method) {\n if (length(TS) != 2) {\n stop('Time series not correct, should be two columns, Date and value.')\n } else if (!grepl('-|/', TS[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.')\n } \n \n \n if (method == 'mon2day') {\n \n data <- apply(TS, MARGIN = 1 , FUN = mon2day)\n \n output <- do.call('rbind', data)\n } else if (method == 'day2mon') {\n Date <- as.Date(TS[, 1])\n year <- format(Date, format = '%Y')\n mon <- format(Date, format = '%m')\n \n data <- aggregate(TS, by = list(mon, year), FUN = mean, na.rm = TRUE)[, 3:4]\n rownames(data) <- 1:dim(data)[1]\n output <- data\n } else {\n stop('method is not correct, check method argument.')\n }\n \n return (output)\n}\n\n#' @importFrom stats aggregate\nresample.list <- function(hyfo, method) {\n checkHyfo(hyfo)\n hyfoData <- hyfo$Data\n Date <- as.POSIXlt(hyfo$Dates$start)\n year <- Date$year + 1900\n mon <- Date$mon + 1\n # hyfoDim <- attributes(hyfoData)$dimensions\n # resample focuses on time dimension. No matter whether the member dimension exists.\n timeIndex <- match('time', attributes(hyfoData)$dimensions)\n dimArray <- 1:length(attributes(hyfoData)$dimensions)\n \n if (method == 'day2mon') {\n hyfoData <- apply(hyfoData, MARGIN = dimArray[-timeIndex], \n function(x) aggregate(x, by = list(mon, year), FUN = mean, na.rm = TRUE)[, 3])\n Date <- aggregate(Date, by = list(mon, year), FUN = mean, na.rm = TRUE)[, 3]\n } else if (method == 'mon2day') {\n message('Under development.')\n }\n \n hyfo$Dates$start <- Date\n hyfo$Data <- hyfoData\n return(hyfo)\n}\n\n\n\n\n#' @importFrom utils tail\n#' @references \n#' \n#' \\itemize{\n#' \\item R Core Team (2015). R: A language and environment for statistical computing. R Foundation for\n#' Statistical Computing, Vienna, Austria. URL http://www.R-project.org/.\n#' }\n#' \nmon2day <- function(monData) {\n Date <- as.Date(monData[1])\n data <- monData[2]\n \n DateY <- format(Date, format = '%Y')\n DateM <- format(Date, format = '%m')\n DateL <- seq(Date, length = 2, by = '1 months')[2] - Date\n \n DateD <- 1:DateL\n \n start <- as.Date(paste(DateY, DateM, DateD[1], sep = '-'))\n end <- as.Date(paste(DateY, DateM, tail(DateD, 1), sep = '-'))\n \n Date <- seq(start, end, by = '1 day')\n \n dailyData <- data.frame(Date = Date, value = rep(data, DateL))\n \n return(dailyData)\n}", - "created" : 1446514556832.000, - "dirty" : false, - "encoding" : "ASCII", - "folds" : "", - "hash" : "2191928295", - "id" : "5ECF1EB5", - "lastKnownWriteTime" : 1446479591, - "path" : "E:/1/R/hyfo/R/resample(generic).R", - "project_path" : "R/resample(generic).R", - "properties" : { - "tempName" : "Untitled1" - }, - "relative_order" : 6, - "source_on_save" : false, - "type" : "r_source" -} \ No newline at end of file diff --git a/.Rproj.user/D53FD3E6/sdb/per/t/C8E82441 b/.Rproj.user/D53FD3E6/sdb/per/t/C8E82441 deleted file mode 100644 index 31d7b34..0000000 --- a/.Rproj.user/D53FD3E6/sdb/per/t/C8E82441 +++ /dev/null @@ -1,18 +0,0 @@ -{ - "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#' @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#' \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#' \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://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 }, S. Herrera \\email{sixto@@predictia.es }\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#' @describeIn biasCorrect\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#' @describeIn 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 <- match('member', attributes(frcData)$dimensions)\n \n # For dataset that has a member part \n if (!is.na(memberIndex)) {\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" : 1446554845419.000, - "dirty" : false, - "encoding" : "ASCII", - "folds" : "", - "hash" : "3817473432", - "id" : "C8E82441", - "lastKnownWriteTime" : 1446431675, - "path" : "E:/1/R/hyfo/R/biasCorrect(generic).R", - "project_path" : "R/biasCorrect(generic).R", - "properties" : { - "tempName" : "Untitled1" - }, - "relative_order" : 7, - "source_on_save" : false, - "type" : "r_source" -} \ No newline at end of file diff --git a/.Rproj.user/D53FD3E6/sdb/per/t/DE19C258 b/.Rproj.user/D53FD3E6/sdb/per/t/DE19C258 index 2906188..fd1fb94 100644 --- a/.Rproj.user/D53FD3E6/sdb/per/t/DE19C258 +++ b/.Rproj.user/D53FD3E6/sdb/per/t/DE19C258 @@ -1,12 +1,12 @@ { - "contents" : "## For package updates information\n\n#' @importFrom utils packageDescription\nhyfoUpdates <- function(){\n page <- readLines('http://yuanchao-xu.github.io/hyfo/')\n updatesLine <- grep('id=\\\\\"updates\"', page)\n versionLine <- updatesLine + 2\n \n version <- unlist(strsplit(page[versionLine], split = ' '))[2]\n version_local <- packageDescription(\"hyfo\")$Version\n \n \n # the first tow digit is the most important part of the version\n version12 <- unlist(strsplit(version, split = \"[.]\"))[1:2]\n version_local12 <- unlist(strsplit(version_local, split = \"[.]\"))[1:2]\n \n sameVersion <- version12 == version_local12\n \n # generate message\n version_msg <- strsplit(strsplit(page[versionLine], split = '

')[[1]][2], split = '

')[[1]]\n infoLine <- versionLine + 2\n info_msg <- strsplit(strsplit(page[infoLine], split = '

')[[1]][2], split = '

')[[1]]\n install_msg <- 'You can update by type in: devtools::install_gihub(\"Yuanchao-Xu/hyfo\")'\n \n message_out <- NULL\n if (any(sameVersion == FALSE)) {\n message_out <- paste(version_msg, info_msg, install_msg, sep = '\\n')\n }\n return(message_out)\n}\n\n.onAttach <- function(libname, pkgname) {\n message_out <- suppressWarnings(try(hyfoUpdates(), silent = TRUE))\n if (!is.null(message_out)) {\n if (grepl('Version', message_out)) {\n packageStartupMessage(message_out)\n }\n }\n}\n", + "contents" : "## For package updates information\n\n#' @importFrom utils packageDescription\nhyfoUpdates <- function(){\n page <- readLines('http://yuanchao-xu.github.io/hyfo/')\n updatesLine <- grep('id=\\\\\"updates\"', page)\n versionLine <- updatesLine + 2\n \n version <- unlist(strsplit(page[versionLine], split = ' '))[2]\n version_local <- packageDescription(\"hyfo\")$Version\n \n \n # the first tow digit is the most important part of the version\n version12 <- unlist(strsplit(version, split = \"[.]\"))[1:2]\n version_local12 <- unlist(strsplit(version_local, split = \"[.]\"))[1:2]\n \n sameVersion <- version12 == version_local12\n \n if (any(sameVersion == FALSE)) {\n # generate message\n version_msg <- strsplit(strsplit(page[versionLine], split = '

')[[1]][2], split = '

')[[1]]\n infoLine_start <- versionLine + 2\n infoLine_end <- grep('

For historical releases and the introduction of updates about each version', page) - 1\n info_msg <- character()\n for (infoLine in infoLine_start:infoLine_end) {\n info_line <- strsplit(strsplit(page[infoLine], split = '>')[[1]][2], split = '<')[[1]][1]\n if (!is.na(info_line)) info_msg <- c(info_msg, info_line)\n }\n \n install_msg <- 'You can update by type in: devtools::install_gihub(\"Yuanchao-Xu/hyfo\")'\n \n message_out <- paste(version_msg, paste(info_msg, collapse = '\\n'), install_msg, sep = '\\n')\n } else message_out <- NULL\n return(message_out)\n}\n\n.onAttach <- function(libname, pkgname) {\n message_out <- suppressWarnings(try(hyfoUpdates(), silent = TRUE))\n if (!is.null(message_out)) {\n if (grepl('Version', message_out)) {\n packageStartupMessage(message_out)\n }\n }\n}\n", "created" : 1446554978988.000, "dirty" : false, "encoding" : "ASCII", "folds" : "", - "hash" : "3670000151", + "hash" : "2271577856", "id" : "DE19C258", - "lastKnownWriteTime" : 1446555739, + "lastKnownWriteTime" : 1446561852, "path" : "E:/1/R/hyfo/R/startup.R", "project_path" : "R/startup.R", "properties" : { diff --git a/NEWS b/NEWS index 59ea925..27d4494 100644 --- a/NEWS +++ b/NEWS @@ -2,8 +2,7 @@ hyfo 1.3.1 ========== Date: 2015.11.3 -- new generic function biasCorrect, extractPeriod, resample, getAnnual, getPreciBar added, - No need to designate inputtype any more, R will detect automatically. +- new generic function biasCorrect, extractPeriod, resample, getAnnual, getPreciBar added. No need to designate inputtype any more, R will detect automatically. - coordinates conversion function extracted. - new user manual about real time bias correction and resample added. diff --git a/R/startup.R b/R/startup.R index 8f87b26..aaa522b 100644 --- a/R/startup.R +++ b/R/startup.R @@ -16,16 +16,21 @@ hyfoUpdates <- function(){ sameVersion <- version12 == version_local12 - # generate message - version_msg <- strsplit(strsplit(page[versionLine], split = '

')[[1]][2], split = '

')[[1]] - infoLine <- versionLine + 2 - info_msg <- strsplit(strsplit(page[infoLine], split = '

')[[1]][2], split = '

')[[1]] - install_msg <- 'You can update by type in: devtools::install_gihub("Yuanchao-Xu/hyfo")' - - message_out <- NULL if (any(sameVersion == FALSE)) { - message_out <- paste(version_msg, info_msg, install_msg, sep = '\n') - } + # generate message + version_msg <- strsplit(strsplit(page[versionLine], split = '

')[[1]][2], split = '

')[[1]] + infoLine_start <- versionLine + 2 + infoLine_end <- grep('

For historical releases and the introduction of updates about each version', page) - 1 + info_msg <- character() + for (infoLine in infoLine_start:infoLine_end) { + info_line <- strsplit(strsplit(page[infoLine], split = '>')[[1]][2], split = '<')[[1]][1] + if (!is.na(info_line)) info_msg <- c(info_msg, info_line) + } + + install_msg <- 'You can update by type in: devtools::install_gihub("Yuanchao-Xu/hyfo")' + + message_out <- paste(version_msg, paste(info_msg, collapse = '\n'), install_msg, sep = '\n') + } else message_out <- NULL return(message_out) }