From a4a0b06eb166e7d9de729cd942f27eeae96bf1f1 Mon Sep 17 00:00:00 2001 From: Vincenzo Coia Date: Mon, 7 Oct 2024 14:34:47 -0700 Subject: [PATCH 1/3] Make a more flexible version of `ch_regime_plot()`. --- DESCRIPTION | 4 +- NAMESPACE | 1 + R/ch_regime_plot2.R | 173 +++++++++++++++++++++++++++++++++++++++++ man/ch_regime_plot2.Rd | 71 +++++++++++++++++ 4 files changed, 248 insertions(+), 1 deletion(-) create mode 100644 R/ch_regime_plot2.R create mode 100644 man/ch_regime_plot2.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 463ae0a..ae23047 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -59,7 +59,9 @@ Imports: tidyhydat, whitebox, datasets, - circular + circular, + vctrs, + rlang Suggests: knitr, testthat, diff --git a/NAMESPACE b/NAMESPACE index b26deff..ab40227 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -43,6 +43,7 @@ export(ch_read_AHCCD_daily) export(ch_read_AHCCD_monthly) export(ch_read_ECDE_flows) export(ch_regime_plot) +export(ch_regime_plot2) export(ch_rfa_distseason) export(ch_rfa_extractamax) export(ch_rfa_julianplot) diff --git a/R/ch_regime_plot2.R b/R/ch_regime_plot2.R new file mode 100644 index 0000000..5a229a8 --- /dev/null +++ b/R/ch_regime_plot2.R @@ -0,0 +1,173 @@ +#' Plots the regime of daily streamflows using quantiles +#' +#' @description Produces a regime hydrograph similar to that in the reference. It shows the flow quantiles for each +#' day of the year and the maximum and minimum. Parameters can be set to change colours and set the y-scale +#' to allow plots of same scale to be produced. +#' +#' @param date,flow Vectors of dates and flow data, or names of the columns +#' containing such data in `data`. +#' @param data Optional data frame containing the data. +#' @param id Gauge code used for creating a smart plot title; could also +#' be a column name in `data`. If `NULL` (the default), title creation +#' is bypassed. +#' @param quant quantiles; default is \code{quant = c(0.95,0.9,0.75,0.5,0.25,0.1,0.05)}. +#' Can be changed but the length must be 7 and the 4th value must be 0.5 (median) +#' @param wyear set \code{wyear = 10} for October, \code{water year = 1} for calendar year, can be any month +#' @param colour if \code{TRUE} plot is in colour, if \code{FALSE} plot is grayscale. +#' @param metadata a data frame of metadata to look for gauge code for `id`, +#' defaults to `HYDAT_list` and is not used if `id` isn't. +#' @param ylab Y axis label. +#' @param ... Other arguments to pass to the `plot()` function. Will take +#' precedence over defaults. +#' +#' @return No value is returned; a standard \R graphic is created. +#' @author Paul Whitfield +#' @importFrom graphics par points polygon legend +#' @importFrom stats quantile +#' @export +#' +#' @references MacCulloch, G. and P. H. Whitfield (2012). Towards a Stream Classification System +#' for the Canadian Prairie Provinces. Canadian Water Resources Journal 37: 311-332. +#' +#' @examples +#' # Refer to columns in the `CAN05AA008` dataset: +#' ch_regime_plot2(Date, Flow, data = CAN05AA008, id = ID) +#' +#' # Optionally, can insert vectors. +#' x <- CAN05AA008$Date +#' y <- CAN05AA008$Flow +#' ch_regime_plot2(x, y) +#' +#' # Override `plot()` defaults using ...; for instance, zero-in on the freshet. +#' ch_regime_plot2(x, y, xlim = c(90, 220)) +ch_regime_plot2 <- function( + date, flow, data = NULL, id = NULL, wyear = 1, colour = TRUE, + metadata = HYDAT_list, + ylab = expression(paste("Mean Daily Discharge (", m^3, "/sec)")), + quant = c(0.95, 0.9, 0.75, 0.5, 0.25, 0.1, 0.05), ...) +{ + flow <- rlang::eval_tidy(rlang::enquo(flow), data = data) + date <- rlang::eval_tidy(rlang::enquo(date), data = data) + v <- vctrs::vec_recycle_common(flow, date) + flow <- v[[1]] + date <- v[[2]] + + id <- unique(rlang::eval_tidy(rlang::enquo(id), data = data)) + + if (is.null(id)) { + title <- NULL + } else { + if (length(id) > 1) { + stop("Received more than one station ID: ", paste0(id, collapse = ", ")) + } + if (length(id) == 0) { + stop("`id` did not evaluate to any gauge code.") + } + sname <- ch_get_wscstation(id, metadata) + title <- sname$Station_lname + } + + + ############################################################################# labels + doy_vals <- ch_doys(date, water_yr = wyear) + year <- doy_vals$year + doy <- doy_vals$doy + + if (wyear != 1) doy <- doy_vals$dwy + + doys <- 366 + doy1 <- c(1:doys) + years <- unique(year) + nyears <- max(years) - min(years) + 1 + min_year <- min(years) - 1 + + ############################################################################# arrays + q <- array(NA, dim = c(nyears, doys)) + + + colr <- c("gray70", "gray50", "gray30", "black", "gray10") + if (colour == TRUE) colr <- c("gray", "cyan", "deepskyblue2", "red", "darkblue") + + ########################################################################## create table of year of daily discharge + for (k in 1:length(year)) { + q[(year[k] - min_year), doy[k]] <- flow[k] + } + + + qquantiles <- quant + qquantiles <- rev(qquantiles) + + regime <- array(NA, dim = c(9, doys)) + + for (jj in 1:doys) { + regime[1, jj] <- min(q[, jj], na.rm = TRUE) + regime[9, jj] <- max(q[, jj], na.rm = TRUE) + for (j in 2:8) { + regime[j, jj] <- stats::quantile(q[, jj], probs = qquantiles[j - 1], na.rm = TRUE) + } + } + + ############################ need to replace Inf and -Inf with NA Infs come from all days being NA + regime[is.infinite(regime)] <- NA + + ########################### create polygons for 0.95-0.05, 0.90-0.1. 0.75-0.25 + mdays <- c(doy1, rev(doy1)) + poly1 <- c(regime[2, ], rev(regime[8, ])) + poly2 <- c(regime[3, ], rev(regime[7, ])) + poly3 <- c(regime[4, ], rev(regime[6, ])) + + ######################################################################### plot start + tscale <- 1.2 + if (!is.null(title) && nchar(title) >= 45) tscale <- 1.0 + if (!is.null(title) && nchar(title) >= 50) tscale <- 0.8 + + # capture plotting parameters, restore on exit + oldpar <- par(no.readonly = TRUE) + on.exit(par(oldpar)) + + par(las = 1) + par(mar = c(3,5,3,1)) + + # Allow override of default plot options by evaluating in a data mask. + plot_defaults <- rlang::env( + xlab = "", xaxt = "n", col = colr[4], + cex = 0.5, ylab = ylab, xlim = c(1, 366), + main = title, cex.main = tscale, + ylim = c(0, max(flow, na.rm = TRUE)) + ) + dots <- rlang::list2(...) + all_args <- union(names(plot_defaults), names(dots)) + names(all_args) <- all_args + all_args <- rlang::parse_exprs(all_args) + all_args_eval <- lapply(all_args, \(x) { + rlang::eval_tidy(x, data = dots, env = plot_defaults) + }) + + # Plot with reconciled arguments + rlang::exec("plot", doy1, regime[9, ], type = "p", !!!all_args_eval) + ch_axis_doy(wyear) + polygon(mdays, poly1, col = colr[1], border = colr[1]) + polygon(mdays, poly2, col = colr[2], border = colr[2]) + polygon(mdays, poly3, col = colr[3], border = colr[3]) + points(doy1, regime[1, ], type = "p", col = colr[4], cex = 0.5) + points(doy1, regime[5, ], type = "l", col = colr[5], lwd = 3) + + ltext1 <- c( + "min / max", + paste( + format(quant[7], nsmall = 2), "-", format(quant[1], nsmall = 2), sep = "" + ), + paste( + format(quant[6], nsmall = 2), "-", format(quant[2], nsmall = 2), sep = "" + ), + paste( + format(quant[5], nsmall = 2), "-", format(quant[3], nsmall = 2), sep = "" + ), + "median" + ) + + lcol1 <- c(colr[4], colr[1], colr[2], colr[3], colr[5]) + legend("topleft", legend = ltext1, col = lcol1, lty = 1, lwd = 3, bty = "n") + ######################################################################### plot end + invisible() +} diff --git a/man/ch_regime_plot2.Rd b/man/ch_regime_plot2.Rd new file mode 100644 index 0000000..e6f6a7d --- /dev/null +++ b/man/ch_regime_plot2.Rd @@ -0,0 +1,71 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ch_regime_plot2.R +\name{ch_regime_plot2} +\alias{ch_regime_plot2} +\title{Plots the regime of daily streamflows using quantiles} +\usage{ +ch_regime_plot2( + date, + flow, + data = NULL, + id = NULL, + wyear = 1, + colour = TRUE, + metadata = HYDAT_list, + ylab = expression(paste("Mean Daily Discharge (", m^3, "/sec)")), + quant = c(0.95, 0.9, 0.75, 0.5, 0.25, 0.1, 0.05), + ... +) +} +\arguments{ +\item{date, flow}{Vectors of dates and flow data, or names of the columns +containing such data in `data`.} + +\item{data}{Optional data frame containing the data.} + +\item{id}{Gauge code used for creating a smart plot title; could also +be a column name in `data`. If `NULL` (the default), title creation +is bypassed.} + +\item{wyear}{set \code{wyear = 10} for October, \code{water year = 1} for calendar year, can be any month} + +\item{colour}{if \code{TRUE} plot is in colour, if \code{FALSE} plot is grayscale.} + +\item{metadata}{a data frame of metadata to look for gauge code for `id`, +defaults to `HYDAT_list` and is not used if `id` isn't.} + +\item{ylab}{Y axis label.} + +\item{quant}{quantiles; default is \code{quant = c(0.95,0.9,0.75,0.5,0.25,0.1,0.05)}. +Can be changed but the length must be 7 and the 4th value must be 0.5 (median)} + +\item{...}{Other arguments to pass to the `plot()` function. Will take +precedence over defaults.} +} +\value{ +No value is returned; a standard \R graphic is created. +} +\description{ +Produces a regime hydrograph similar to that in the reference. It shows the flow quantiles for each +day of the year and the maximum and minimum. Parameters can be set to change colours and set the y-scale +to allow plots of same scale to be produced. +} +\examples{ +# Refer to columns in the `CAN05AA008` dataset: +ch_regime_plot2(Date, Flow, data = CAN05AA008, id = ID) + +# Optionally, can insert vectors. +x <- CAN05AA008$Date +y <- CAN05AA008$Flow +ch_regime_plot2(x, y) + +# Override `plot()` defaults using ...; for instance, zero-in on the freshet. +ch_regime_plot2(x, y, xlim = c(90, 220)) +} +\references{ +MacCulloch, G. and P. H. Whitfield (2012). Towards a Stream Classification System +for the Canadian Prairie Provinces. Canadian Water Resources Journal 37: 311-332. +} +\author{ +Paul Whitfield +} From 4c61f6ccabfb68feb5601a1e6b941dcff291c23c Mon Sep 17 00:00:00 2001 From: Vincenzo Coia Date: Wed, 16 Oct 2024 14:30:21 -0700 Subject: [PATCH 2/3] add scale_x/y_gumbel --- R/scale_gumbel.R | 67 ++++++++++++++++++++++++++++++++++++++++++ man/gumbelAEP_trans.Rd | 16 ++++++++++ man/gumbelRP_trans.Rd | 16 ++++++++++ man/gumbel_spacing.Rd | 45 ++++++++++++++++++++++++++++ 4 files changed, 144 insertions(+) create mode 100644 R/scale_gumbel.R create mode 100644 man/gumbelAEP_trans.Rd create mode 100644 man/gumbelRP_trans.Rd create mode 100644 man/gumbel_spacing.Rd diff --git a/R/scale_gumbel.R b/R/scale_gumbel.R new file mode 100644 index 0000000..f68e457 --- /dev/null +++ b/R/scale_gumbel.R @@ -0,0 +1,67 @@ +#' Gumbel-transformed axes +#' +#' Transforms return period or annual exceedance probability to be spaced +#' according to a reduced Gumbel variate. Can be added as a scale layer to +#' a ggplot2 graphic. +#' +#' @param ... Arguments to pass to `scale_x_continuous` or `scale_y_continuous` +#' from ggplot2. +#' @return The same output as `scale_x_continuous` or `scale_y_continuous`, +#' but with the appropriate Gumbel spacing. +#' @rdname gumbel_spacing +#' @examples +#' library(ggplot2) +#' df <- data.frame( +#' y = sort(stats::rexp(100), decreasing = TRUE), +#' aep = 1:100 / 101 +#' ) +#' +#' ggplot(df, aes(aep, y)) + +#' geom_point() + +#' scale_x_gumbelAEP() +#' +#' ggplot(df, aes(1 / aep, y)) + +#' geom_point() + +#' scale_x_gumbelRP("Return Period") +#' @export +scale_x_gumbelRP <- function(...) { + ggplot2::scale_x_continuous(..., trans = gumbelRP_trans) +} + +#' @rdname gumbel_spacing +#' @export +scale_y_gumbelRP <- function(...) { + ggplot2::scale_y_continuous(..., trans = gumbelRP_trans) +} + +#' @rdname gumbel_spacing +#' @export +scale_x_gumbelAEP <- function(...) { + ggplot2::scale_x_continuous(..., trans = gumbelAEP_trans) +} + +#' @rdname gumbel_spacing +#' @export +scale_y_gumbelAEP <- function(...) { + ggplot2::scale_y_continuous(..., trans = gumbelAEP_trans) +} + +#' Gumbel transformations used for ggplot2 scales +#' @export +gumbelRP_trans <- scales::trans_new( + "gumbelRP", + transform = function(x) -log(-log(1 - 1 / x)), + inverse = function(x) 1 / (1 - exp(-exp(-x))), + breaks = scales::breaks_log(), + domain = c(1, Inf) +) + +#' Gumbel transformations used for ggplot2 scales +#' @export +gumbelAEP_trans <- scales::trans_new( + "gumbelAEP", + transform = function(x) -log(-log(1 - x)), + inverse = function(x) 1 - exp(-exp(-x)), + breaks = scales::breaks_log(), + domain = c(1e-100, 1) +) \ No newline at end of file diff --git a/man/gumbelAEP_trans.Rd b/man/gumbelAEP_trans.Rd new file mode 100644 index 0000000..a71fe4b --- /dev/null +++ b/man/gumbelAEP_trans.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/scale_gumbel.R +\docType{data} +\name{gumbelAEP_trans} +\alias{gumbelAEP_trans} +\title{Gumbel transformations used for ggplot2 scales} +\format{ +An object of class \code{transform} of length 9. +} +\usage{ +gumbelAEP_trans +} +\description{ +Gumbel transformations used for ggplot2 scales +} +\keyword{datasets} diff --git a/man/gumbelRP_trans.Rd b/man/gumbelRP_trans.Rd new file mode 100644 index 0000000..54b1307 --- /dev/null +++ b/man/gumbelRP_trans.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/scale_gumbel.R +\docType{data} +\name{gumbelRP_trans} +\alias{gumbelRP_trans} +\title{Gumbel transformations used for ggplot2 scales} +\format{ +An object of class \code{transform} of length 9. +} +\usage{ +gumbelRP_trans +} +\description{ +Gumbel transformations used for ggplot2 scales +} +\keyword{datasets} diff --git a/man/gumbel_spacing.Rd b/man/gumbel_spacing.Rd new file mode 100644 index 0000000..06b2b98 --- /dev/null +++ b/man/gumbel_spacing.Rd @@ -0,0 +1,45 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/scale_gumbel.R +\name{scale_x_gumbelRP} +\alias{scale_x_gumbelRP} +\alias{scale_y_gumbelRP} +\alias{scale_x_gumbelAEP} +\alias{scale_y_gumbelAEP} +\title{Gumbel-transformed axes} +\usage{ +scale_x_gumbelRP(...) + +scale_y_gumbelRP(...) + +scale_x_gumbelAEP(...) + +scale_y_gumbelAEP(...) +} +\arguments{ +\item{...}{Arguments to pass to `scale_x_continuous` or `scale_y_continuous` +from ggplot2.} +} +\value{ +The same output as `scale_x_continuous` or `scale_y_continuous`, +but with the appropriate Gumbel spacing. +} +\description{ +Transforms return period or annual exceedance probability to be spaced +according to a reduced Gumbel variate. Can be added as a scale layer to +a ggplot2 graphic. +} +\examples{ +library(ggplot2) +df <- data.frame( + y = sort(stats::rexp(100), decreasing = TRUE), + aep = 1:100 / 101 +) + +ggplot(df, aes(aep, y)) + + geom_point() + + scale_x_gumbelAEP() + +ggplot(df, aes(1 / aep, y)) + + geom_point() + + scale_x_gumbelRP("Return Period") +} From d66d474e6754289d29d9417bf2e8b592403e7048 Mon Sep 17 00:00:00 2001 From: Vincenzo Coia Date: Wed, 16 Oct 2024 14:31:10 -0700 Subject: [PATCH 3/3] make ch_ams_timeseries() --- NAMESPACE | 7 ++++++ R/ch_ams_timeseries.R | 54 ++++++++++++++++++++++++++++++++++++++++ man/ch_ams_timeseries.Rd | 44 ++++++++++++++++++++++++++++++++ 3 files changed, 105 insertions(+) create mode 100644 R/ch_ams_timeseries.R create mode 100644 man/ch_ams_timeseries.Rd diff --git a/NAMESPACE b/NAMESPACE index ab40227..e50b516 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -9,6 +9,7 @@ S3method(ch_rfa_extractamax,formula) S3method(ch_rfa_seasonstat,data.frame) S3method(ch_rfa_seasonstat,default) S3method(ch_rfa_seasonstat,formula) +export(ch_ams_timeseries) export(ch_axis_doy) export(ch_binned_MannWhitney) export(ch_booth_plot) @@ -69,6 +70,12 @@ export(ch_wbt_flow_direction) export(ch_wbt_pourpoints) export(ch_wbt_removesinks) export(ch_wtr_yr) +export(gumbelAEP_trans) +export(gumbelRP_trans) +export(scale_x_gumbelAEP) +export(scale_x_gumbelRP) +export(scale_y_gumbelAEP) +export(scale_y_gumbelRP) import(circular) importFrom(Kendall,MannKendall) importFrom(dplyr,left_join) diff --git a/R/ch_ams_timeseries.R b/R/ch_ams_timeseries.R new file mode 100644 index 0000000..7b4768a --- /dev/null +++ b/R/ch_ams_timeseries.R @@ -0,0 +1,54 @@ +#' Produce a Time Series Plot +#' +#' Uses ggplot2 to produce a time series plot of values (e.g., flow) +#' over time, where time is binned by year. Subsequent years are connected with +#' a line, which is broken when data are missing. +#' +#' @param date Vector of dates or years, or name of data column containing the +#' dates (unquoted). +#' @param flow Vector of flows, or name of data column containing the flows +#' (unquoted). +#' @param data Optional; data frame containing the specified data. If +#' specified, this data frame will be searched first when looking for data +#' vectors. +#' @returns A ggplot2 object. This means that you are able to add more layers +#' downstream. +#' @examples +#' library(lubridate) +#' set.seed(42) +#' dates <- ymd(paste( +#' c(1991:1993, 1995, 1997:2012), "-", +#' sample(1:12, size = 20, replace = TRUE), "-", +#' sample(1:28, size = 20, replace = TRUE) +#' )) +#' y <- stats::rexp(20) +#' y[18] <- NA +#' ch_ams_timeseries(dates, y) +#' df <- data.frame(year = year(dates), flow = y) +#' ch_ams_timeseries(year, 35.31467 * flow, data = df) + +#' ylab("Peak Instantaneous Flow (cfs)") + +#' theme_bw() +#' @export +ch_ams_timeseries <- function(date, flow, data = NULL) { + in_date <- rlang::enquo(date) + in_flow <- rlang::enquo(flow) + name_flow <- names(rlang::quos_auto_name(list(in_flow))) + date <- rlang::eval_tidy(in_date, data = data) + flow <- rlang::eval_tidy(in_flow, data = data) + if (lubridate::is.Date(date) || lubridate::is.POSIXt(date)) { + year <- lubridate::year(date) + } else { + year <- date + } + xy <- vctrs::vec_recycle_common(year, flow) # Enforce strict recycling. + year <- xy[[1]] + flow <- xy[[2]] + complete_years <- seq(min(year, na.rm = TRUE), max(year, na.rm = TRUE)) + complete_df <- data.frame(year = complete_years) + original_df <- data.frame(year = year, flow = flow) + df <- dplyr::left_join(complete_df, original_df, by = "year") + if (nrow(df) == 0) return(ggplot2::ggplot()) + ggplot2::ggplot(df, ggplot2::aes(year, flow)) + + ggplot2::geom_point() + + ggplot2::geom_line() +} diff --git a/man/ch_ams_timeseries.Rd b/man/ch_ams_timeseries.Rd new file mode 100644 index 0000000..8b25b9e --- /dev/null +++ b/man/ch_ams_timeseries.Rd @@ -0,0 +1,44 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ch_ams_timeseries.R +\name{ch_ams_timeseries} +\alias{ch_ams_timeseries} +\title{Produce a Time Series Plot} +\usage{ +ch_ams_timeseries(date, flow, data = NULL) +} +\arguments{ +\item{date}{Vector of dates or years, or name of data column containing the +dates (unquoted).} + +\item{flow}{Vector of flows, or name of data column containing the flows +(unquoted).} + +\item{data}{Optional; data frame containing the specified data. If +specified, this data frame will be searched first when looking for data +vectors.} +} +\value{ +A ggplot2 object. This means that you are able to add more layers +downstream. +} +\description{ +Uses ggplot2 to produce a time series plot of values (e.g., flow) +over time, where time is binned by year. Subsequent years are connected with +a line, which is broken when data are missing. +} +\examples{ +library(lubridate) +set.seed(42) +dates <- ymd(paste( + c(1991:1993, 1995, 1997:2012), "-", + sample(1:12, size = 20, replace = TRUE), "-", + sample(1:28, size = 20, replace = TRUE) +)) +y <- stats::rexp(20) +y[18] <- NA +ch_ams_timeseries(dates, y) +df <- data.frame(year = year(dates), flow = y) +ch_ams_timeseries(year, 35.31467 * flow, data = df) + + ylab("Peak Instantaneous Flow (cfs)") + + theme_bw() +}