Nothing
#' Change point time series plot
#'
#' Compiles change point information for all metrics and outputs a daily flow
#' time series plot overlain with a bar plot of changepoint counts by year.
#' @param metrics output from \code{\link{metrics.all}}
#' @param type character indicating which type of metric to compile change points for.
#' Options are "h" for high flow metrics, "l" for low flow metrics, "b" for baseflow
#' metrics, or "a" for all 30 metrics (10 high, 10 low, 10 baseflow).
#' @param text optional character string for margin text, e.g. for station name,
#' location, or other notes. Set to NULL if no margin text is wanted, or set to "d"
#' to use default text containing the station ID, station name, and province/state
#' returned from \code{\link{station.info}}.
#' @return When type="a", returns a data.frame of changepoint counts by metric
#' type and year.
#' @author Jennifer Dierauer
#' @seealso \code{\link{metrics.all}}
#' @export
#' @examples
#' # load results from metrics.all function for the Caniapiscau River
#' data(caniapiscau.res)
#'
#' # plot changepoints for all metrics
#' screen.cpts(caniapiscau.res, type="l")
screen.cpts <- function(metrics, type="a", text=NULL) {
opar <- graphics::par(no.readonly = TRUE)
res <- metrics[[2]]
TS <- metrics[[3]]
cptsh <- list()
cptsl <- list()
cptsb <- list()
cmcolors <- grDevices::cm.colors(20, alpha=0.5)
colh <- cmcolors[1]
coll <- cmcolors[20]
colb <- cmcolors[10]
if (type == "h") {
mtitle <- "Changepoints in High Flow Metrics"
rylab <- "Number of Changepoints (max = 10)"
types <- "Change Point Count"
mcol <- colh
}
if (type == "l") {
mtitle <- "Changepoints in Low Flow Metrics"
rylab <- "Number of Changepoints (max = 10)"
types <- "Change Point Count"
mcol <- coll
}
if (type == "b") {
mtitle <- "Changepoints in Baseflow Metrics"
rylab <- "Number of Changepoints (max = 10)"
types <- "Change Point Count"
mcol <- colb
}
if (type == "a") {
mtitle <- "Changepoints in All Metrics"
rylab <- "Number of Changepoints (max = 30)"
types <- c("Baseflow Changepoints", "Low Flow Changepoints",
"High Flow Changepoints")
mcol <- c(colb, coll, colh)
}
# put all change points in a list
for (i in 1:10) {
metname <- names(res[i])
cptsh[[metname]] <- res[[i]]$cpts
}
for (i in 11:20) {
metname <- names(res[i])
cptsl[[metname]] <- res[[i]]$cpts
}
for (i in 21:30) {
metname <- names(res[i])
cptsb[[metname]] <- res[[i]]$cpts
}
## identify start and end dates for plot, set axis limits
Year1 <- TS$year[1]
Start <- as.Date(paste(Year1, "-01-01", sep=""))
YearLast <- TS$year[length(TS$year)]
End <- as.Date(paste(YearLast, "-12-31", sep=""))
myxlims <- c(Start, End)
myylims <- c(0, 1.1 * max(TS$Flow))
NumYears <- as.numeric(YearLast) - as.numeric(Year1) + 1
Year_List <- seq(from=Year1, to=YearLast, by=1)
CountsH <- c(numeric(0))
CountsL <- c(numeric(0))
CountsB <- c(numeric(0))
for (i in 1:length(cptsh)) {
cpts.sub <- cptsh[[i]]
if (!is.na(cpts.sub[1])) {
cpts.sub <- attr(cpts.sub, "times")
for (j in 1:length(cpts.sub)) {
ind <- as.numeric(substr(as.character(cpts.sub[j]), 1, 4))
CountsH <- c(CountsH, ind)
}
}
}
for (i in 1:length(cptsl)) {
cpts.sub <- cptsl[[i]]
if (!is.na(cpts.sub[1])) {
cpts.sub <- attr(cpts.sub, "times")
for (j in 1:length(cpts.sub)) {
ind <- as.numeric(substr(as.character(cpts.sub[j]), 1, 4))
CountsL <- c(CountsL, ind)
}
}
}
for (i in 1:length(cptsb)) {
cpts.sub <- cptsb[[i]]
if (!is.na(cpts.sub[1])) {
cpts.sub <- attr(cpts.sub, "times")
for (j in 1:length(cpts.sub)) {
ind <- as.numeric(substr(as.character(cpts.sub[j]), 1, 4))
CountsB <- c(CountsB, ind)
}
}
}
y1 <- expression(paste("Discharge (m" ^{3}, "/s)"))
graphics::par(mar = c(3,5,2,5))
if (!is.null(text)) {graphics::par(oma=c(0,0,1,0))}
graphics::plot(TS$Date, TS$Flow, pch=19, cex=0.3, ylab=y1, xlab="",
xlim=myxlims, ylim=myylims, col="grey50")
graphics::title(main=mtitle)
if (!is.null(text)) {
if (text == "d") {
stinfo <- station.info(TS, Plot=F)
text <- paste("ID: ", stinfo[1],", NAME: ", stinfo[2],
", PROV/STATE: ", stinfo[3], sep = "")
}
graphics::mtext(text, side=3, line=0, outer=T, cex=0.7)
}
### add shaded polygons covering missing data periods
MissingDays <- NA.runs(TS)
polycol <- grDevices::gray(0.5, alpha=0.5)
if (length(MissingDays$Start) != 0) {
for (i in 1:length(MissingDays$Start)) {
xx <- c(MissingDays$Start[i], MissingDays$Start[i],
MissingDays$End[i], MissingDays$End[i])
yy <- c(-100, 1.5*max(TS$Flow), 1.5*max(TS$Flow), -100)
graphics::polygon(xx,yy,col=polycol, border=NA)
}
}
CountsH <- CountsH[CountsH < YearLast - 1]
CountsB <- CountsB[CountsB < YearLast - 1]
CountsL <- CountsL[CountsL < YearLast - 1]
hv1 <- rep(0, length(Year_List))
hv2 <- rep(0, length(Year_List))
hv3 <- rep(0, length(Year_List))
for (i in 1:length(Year_List)) {
hv1[i] <- length(CountsH[CountsH == Year_List[i]])
hv2[i] <- length(CountsL[CountsL == Year_List[i]])
hv3[i] <- length(CountsB[CountsB == Year_List[i]])
}
## add histogram of changepoints
if (type == "a") {
graphics::par(new = T)
graphics::barplot(rbind(hv3, hv2, hv1), col=mcol, xlab="", ylim=c(0, 30),
xaxt="n", yaxt="n")
}
if (type == "h") {
graphics::par(new = T)
graphics::barplot(hv1, col=mcol, xlab="", ylim=c(0, 10),
xaxt="n", yaxt="n")
}
if (type == "l") {
graphics::par(new = T)
graphics::barplot(hv2, col=mcol, xlab="", ylim=c(0, 10),
xaxt="n", yaxt="n")
}
if (type == "b") {
graphics::par(new = T)
graphics::barplot(hv3, col=mcol, xlab="", ylim=c(0, 10),
xaxt="n", yaxt="n")
}
graphics::axis(side=4)
graphics::mtext(side = 4, line=3, rylab)
graphics::legend("topright", legend = c("Missing Data Periods", types), col=c(polycol, mcol)
, fill=c(polycol, mcol), bty="n")
res <- data.frame(Year=Year_List, HighFlow=hv1, LowFlow=hv2, Baseflow=hv3)
if (type == "a") {return(res)}
# restore settings on exit
on.exit(suppressWarnings(graphics::par(opar)))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.