R/PCP_plot.R

In sharpshootR: A Soil Survey Toolkit

# cumulative PPT within a single water year
# values are pre-sorted
.cumulativePPT <- function(i) {
  
  # replace NA with 0
  i$value <- ifelse(is.na(i$value), 0, i$value)
  
  ## TODO: think about how to represent flashy-ness
  # use first derivative
  # i$value <- c(0, diff(i$value))
  
  # cumulative PPT
  i$cumulative_ppt <- cumsum(i$value)
  
  # number of days in summation
  i$n <- nrow(i)
  
  ## TODO: subset to rows of interest
  return(i)
}

# compute cumulative PPT by water year
.prepareCumulativePPT <- function(d) {
  ## NOTE: requires sharpshootR >= 1.4.02
  # re-order just in case
  d <- d[order(d$water_year, d$water_day), ]
  
  # compute cumulative PPT by water year, ordered by water day
  dd <- lapply(split(d, d$water_year), .cumulativePPT)
  dd <- do.call('rbind', dd)
  
  return(dd)
}



# get years and total PPT which are closest to given percentiles of annual PPT
.exemplarYears <- function(d) {
  
  # annual PPT by water year
  ppt.by.wy <- tapply(d$cumulative_ppt, d$water_year, max, na.rm = TRUE)
  
  # interesting percentiles of annual PPT
  ppt.q <- quantile(ppt.by.wy, probs = c(0.05, 0.25, 0.5, 0.75, 0.95), na.rm = TRUE)
  
  # find exemplars for select percentiles
  ppt.abs.diff <- abs(outer(ppt.by.wy, ppt.q, FUN = "-"))
  ppt.closest.yrs <- apply(ppt.abs.diff, MARGIN = 2, FUN = which.min)
  
  # this is a named vector of annual PPT (names are water_year)
  exemplar.data <- ppt.by.wy[ppt.closest.yrs]
  
  return(exemplar.data)
}


## TODO: generalize to other sources of data: DAYMET / SCAN / SNOTEL / Henry / etc.
## TODO: add percentiles-by-water-day method
## TODO: this.year should default to the last year in the series
## TODO: explain differences between exemplar years and boxplot at current water day -- not the same

## TODO: remove current WY from exemplars if not at the end of the year

#' @title Percentiles of Cumulative Precipitation
#' 
#' @description Generate a plot representing percentiles of cumulative precipitation, given a historic record, and criteria for selecting a year of data for comparison.
#'
#' @param x result from `CDECquery` for now, will need to generalize to other sources
#' @param this.year a single water year, e.g. 2020
#' @param this.day optional integer representing days since start of selected water year
#' @param method 'exemplar' or 'daily', currently 'exemplar' is the only method available
#' @param q.color color of percentiles cumulative precipitation
#' @param c.color color of selected year
#' @param ... additional arguments to `plot`
#' 
#' @details This is very much a work in progress. Further examples at \url{https://ncss-tech.github.io/AQP/sharpshootR/CDEC.html}, and \url{https://ncss-tech.github.io/AQP/sharpshootR/cumulative-PPT.html}.
#'
#' @return nothing, this function is called to create graphical output
#' @author D.E. Beaudette
#' 
#' @seealso [soilDB::waterDayYear()]
#' @keywords hplots
#' 
#' @export
#'
PCP_plot <- function(x, this.year, this.day = NULL, method = 'exemplar', q.color = 'RoyalBlue', c.color = 'firebrick', ...) {
  
  # hack for R CMD check
  cumulative_ppt <- NULL
  
  # sanity check
  method <- match.arg(method)
  
  # water year range
  wy.range <- range(x$water_year)
  
  # prepare data for plotting / extract summaries
  xx <- .prepareCumulativePPT(x)
  e <- .exemplarYears(xx)
  
  # convenience objects for plotting 
  exemplar.yrs <- as.numeric(names(e))
  
  if(is.null(this.day)) {
    # subset year
    idx <- which(xx$water_year == this.year)
    this.year.data <- xx[idx, ]  
  } else {
    # truncate to interval {2, 365}
    this.day <- pmax(this.day, 2)
    this.day <- pmin(this.day, 365)
    # subset year and days <= selected water day
    idx <- which(xx$water_year == this.year & xx$water_day <= this.day)
    this.year.data <- xx[idx, ]
  }
  
  
  ## TODO: warn / remove years with less than 330 (?) days of data
  # tapply(xx$water_day, xx$water_year, length)
  
  
  # current year positional elements
  # last real date in series
  mrd <- as.character(max(this.year.data$datetime, na.rm = TRUE))
  # last water date in series
  mwd <- max(this.year.data$water_day, na.rm = TRUE)
  # last cumulative PPT in series
  mcp <- max(this.year.data$cumulative_ppt, na.rm = TRUE)
  
  # water-day stats with PPT < 0.1
  dry.wd <- xx$water_day[xx$cumulative_ppt < 0.1]
  dry.wd.q <- quantile(dry.wd, probs = c(0.05, 0.25, 0.5, 0.75, 0.95), na.rm = TRUE)
  
  # PPT stats on current water-day
  this.wd.data <- xx$cumulative_ppt[xx$water_day == mwd]
  this.wd.data.q <- quantile(this.wd.data, probs = c(0.05, 0.25, 0.5, 0.75, 0.95), na.rm = TRUE)
  
  # reasonable axis in "real" dates
  date.axis <- data.frame(
    d = seq.Date(from = as.Date('2000-10-01'), to = as.Date('2001-09-30'), by = '2 weeks')
  )
  
  # integrate water year / day
  wyd <- waterDayYear(date.axis$d)
  date.axis$wr <- wyd$wy
  date.axis$wd <- wyd$wd
  
  # generic labeling of months and days
  date.axis$lab <- format(date.axis$d, "%d\n%b")
  
  # horizontal grid lines
  y.grid <- pretty(xx$cumulative_ppt, n = 10)
  
  # keep track of water years, label x-axis with this
  xlab.text <- sprintf("Water Years %s - %s", wy.range[1], wy.range[2])
  
  # all data, establish plot area
  plot(cumulative_ppt ~ water_day, data = xx, col=grey(0.9), type='n', axes=FALSE, xlim=c(-2, 375), ylim=c(-2, max(xx$cumulative_ppt, na.rm = TRUE)), xlab = xlab.text, ... )
  # plot(cumulative_ppt ~ water_day, data=xx, col=grey(0.9), type='n', axes=FALSE, xlim=c(-5, 370), xlab=xlab.text)
  
  # grid
  abline(h = y.grid, v = date.axis$wd, col = 'lightgray', lty = 3)
  
  
  # filter trivial rainfall at the start of the year
  # filter +1 day due to leap year
  idx <- which(xx$cumulative_ppt >= 0.1 & xx$water_day <= 365, )
  vars <- c('water_day', 'cumulative_ppt', 'water_year')
  xx.q <- xx[idx, vars]
  
  # exemplar years (percentiles) form boundary of polygons
  xx.q05 <- xx.q[xx.q$water_year == exemplar.yrs[1], ]
  xx.q25 <- xx.q[xx.q$water_year == exemplar.yrs[2], ]
  xx.q50 <- xx.q[xx.q$water_year == exemplar.yrs[3], ]
  xx.q75 <- xx.q[xx.q$water_year == exemplar.yrs[4], ]
  xx.q95 <- xx.q[xx.q$water_year == exemplar.yrs[5], ]
  
  # q05--q95 polygon
  p.1.x <- c(xx.q05$water_day, rev(xx.q95$water_day))
  p.1.y <- c(xx.q05$cumulative_ppt, rev(xx.q95$cumulative_ppt))
  p.1.col <- rgb(t(col2rgb(q.color)), maxColorValue = 255, alpha = 33)
  
  # q25--q95 polygon
  p.2.x <- c(xx.q25$water_day, rev(xx.q75$water_day))
  p.2.y <- c(xx.q25$cumulative_ppt, rev(xx.q75$cumulative_ppt))
  p.2.col <- rgb(t(col2rgb(q.color)), maxColorValue = 255, alpha = 100)
  
  # transparency is visual cue
  polygon(x = p.1.x, y = p.1.y, col = p.1.col, border = q.color)
  polygon(x = p.2.x, y = p.2.y, col = p.2.col, border = q.color)
  
  # q50
  lines(cumulative_ppt ~ water_day, data=xx[xx$water_year == exemplar.yrs[3], ],
        subset=cumulative_ppt >= 0.1, lwd=2, lty=1, col=q.color, type='l')
  
  # current year, non-zero values only
  lines(cumulative_ppt ~ water_day, data=this.year.data, 
        subset=cumulative_ppt >= 0.1, lwd=2, col=c.color, type='l')
  
  # add axes
  axis(side=1, at = date.axis$wd, labels = date.axis$lab, cex.axis=0.55, las=1)
  axis(side=2, las=1, at = y.grid, cex.axis=0.75)
  
  # annotate exemplar years
  txt <- sprintf("%s (%s)", as.character(round(e)), as.character(exemplar.yrs))
  text(x = 362, y=e, labels = txt, pos=4, cex=0.65, font=2)
  # text(x = 370, y=e+2, labels = , adj = c(0,0), cex=0.65, font=2)
  
  # add current year's cumulative PPT
  points(x = mwd, y = mcp, pch = 22, bg = c.color)
  
  ## TODO: this threshold and position of bxp works for units of inches but not mm
  # add boxplot of starting water-day where cumulative PPT > 0.1
  bxp.data <- list(stats=matrix(dry.wd.q, ncol=1), n=length(dry.wd.q), out=NULL, group=1, names="")
  # add to current plot
  bxp(bxp.data, at=-1, add=TRUE, show.names=FALSE, outline = FALSE, axes = FALSE, boxwex=2, border=q.color, horizontal = TRUE)
  
  # add box plot of cumulative PPT on current water-day
  if(mcp > 0) {
    # add percentiles for current water day via boxplot
    # create data needed by bxp()
    bxp.data <- list(stats=matrix(this.wd.data.q, ncol=1), n=length(this.wd.data), out=NULL, group=1, names="")
    # add to current plot
    bxp(bxp.data, at=0, add=TRUE, show.names=FALSE, outline = FALSE, axes = FALSE, boxwex=10, border=c.color)
    # annotate with customized quantiles
    text(x=0, y=this.wd.data.q, labels = names(this.wd.data.q), cex=0.65, pos=2, col=c.color)
    
    # annotate with current cumulative PPT
    text(x=0, y=mcp, labels = round(mcp, 1), col=c.color, font=2, cex=0.75, pos=4)
  }
  
  
  # helper lines confusing if last water day < ~ 25
  if(mwd > 25) {
    # helper lines confusing if last cumulative PPT < 1
    if(mcp > 0) {
      # helper lines
      segments(x0 = 20, y0 = mcp, x1 = mwd, y1 = mcp, col = alpha(c.color, 0.5))
      segments(x0 = mwd, y0 = mcp, x1 = mwd, y1 = 1, col = alpha(c.color, 0.5))
    }
    
    # annotate current (real) date
    text(x=mwd, y=-2.25, labels = mrd, font=2, cex=0.75, col=c.color)  
  }
  
  
  # re-make inner polygon color for legend
  p.2.col.leg <- rgb(t(col2rgb(q.color)), maxColorValue = 255, alpha = 150)
  
  # legend
  legend(
    'top', 
    legend = c('5th--95th Percentiles', '25th--75th Percentiles', '50th Percentile', sprintf('Current (%s)', this.year)), 
    pt.bg = c(p.1.col, p.2.col.leg, q.color, c.color), 
    pch = 22, 
    pt.cex = 1.25, 
    cex = 0.75, 
    ncol = 4, 
    bg = 'white'
  )
  
  # TODO: return information used to make figure
}


## alternative approach: percentiles by water day, over all water years
## NOT monotonic functions

# # TODO: would be nice to know number of years
# # TODO: filtering on number of days in summation?
# # percentiles over all years, by water day
# .PPT_pctiles <- function(i) {
#   res <- quantile(i$cumulative_ppt, probs=c(0.05, 0.25, 0.5, 0.75, 0.95), na.rm = TRUE)
#   return(res)
# }
# 
# ## not quite right, leap years aren't adequately handled     
# xx <- ddply(x, 'water_year', .fun=.cumulativePPT)
# xxx <- ddply(xx, 'water_day', .fun=.PPT_pctiles)  
# 
# zz <- subset(xx, subset=water_year == 2019)
# zzz <- ddply(zz, 'water_day', .fun=.PPT_pctiles)  
# 
# # almost right
# matplot(xxx[1:365, 1], xxx[1:365, -1], type='l', col='RoyalBlue', 
#         lwd=c(1,1,2,1,1), lty=c(3,2,1,2,3), 
#         las=1, xlab='Water Day', ylab='Cumulative PPT (inches)', 
#         main='Sonora Ranger Station (SOR)\n1981--2019')
# 
# matlines(zzz[1:365, 1], zzz[1:365, 4], lwd=2, col='firebrick')