This document illustrates how to create stormwater event plots showing precipitation and streamflow. These plots are used to evaluate runoff response times and durations during storm events.
library(devtools)
library(lubridate)
library(plyr)
library(ggplot2)
library(gridExtra)
These functions are stored as gists on Github for common usage. Note that this is the same as downloading these files and running source().
source_gist("https://gist.github.com/walkerjeffd/6179773") # zoo.regular()
source_gist("https://gist.github.com/walkerjeffd/6179768") # assign.events() and events.summary()
Hourly precipitation data were obtained from the NRCC and stored in a CSV file.
# load precipitation data
pcp <- read.csv("..//data//precip//logan.csv", as.is = TRUE)
pcp <- transform(pcp, Date = mdy(Date), DateTime = mdy_hm(DateTime)) # parse dates
# convert to regular zoo object (i.e. continuous time series that has no
# missing time steps)
pcp.zoo <- zoo.regular(pcp$DateTime, pcp$Precip)
# convert back to dataframe
pcp.df <- data.frame(DATETIME = time(pcp.zoo), PRECIP = coredata(pcp.zoo))
head(pcp.df)
## DATETIME PRECIP
## 1 2010-12-31 19:00:00 0
## 2 2010-12-31 20:00:00 0
## 3 2010-12-31 21:00:00 0
## 4 2010-12-31 22:00:00 0
## 5 2010-12-31 23:00:00 0
## 6 2011-01-01 00:00:00 0
15-minute streamflow data were downloaded from the USGS and stored in a CSV file.
q <- read.csv("..//data//usgs//01103025_15min_2011_2014.csv", as.is = TRUE)
q$Datetime <- mdy_hm(q$Datetime) # parse datetimes
q$Datetime <- q$Datetime - hours(ifelse(q$TZ == "EDT", 1, 0)) # convert from EST/EDT to EST only
q <- q[order(q$Datetime), ] # sort by datetime
head(q)
## Datetime TZ Stage Velocity Flow
## 1 2011-04-14 23:00:00 EDT 2.28 0.50 17
## 2 2011-04-14 23:15:00 EDT 2.28 0.53 18
## 3 2011-04-14 23:30:00 EDT 2.29 0.53 18
## 4 2011-04-14 23:45:00 EDT 2.30 0.50 17
## 5 2011-04-15 00:00:00 EDT 2.30 0.50 17
## 6 2011-04-15 00:15:00 EDT 2.30 0.39 13
Convert the 15-minute streamflow data to hourly data.
q.hr <- ddply(q, .(Datetime = floor_date(Datetime, unit = "hour")), summarize,
Stage = mean(Stage), Velocity = mean(Velocity), Flow = mean(Flow))
head(q.hr)
## Datetime Stage Velocity Flow
## 1 2011-04-14 23:00:00 2.287 0.5150 17.50
## 2 2011-04-15 00:00:00 2.308 0.4200 14.25
## 3 2011-04-15 01:00:00 2.333 0.4825 16.50
## 4 2011-04-15 02:00:00 2.367 0.4200 14.75
## 5 2011-04-15 03:00:00 2.375 0.4025 14.00
## 6 2011-04-15 04:00:00 2.373 0.4450 15.25
df <- merge(pcp.df, q.hr, by.x = "DATETIME", by.y = "Datetime", all = FALSE)
names(df) <- toupper(names(df))
head(df)
## DATETIME PRECIP STAGE VELOCITY FLOW
## 1 2011-04-14 19:00:00 0 2.287 0.5150 17.50
## 2 2011-04-14 20:00:00 0 2.308 0.4200 14.25
## 3 2011-04-14 21:00:00 0 2.333 0.4825 16.50
## 4 2011-04-14 22:00:00 0 2.367 0.4200 14.75
## 5 2011-04-14 23:00:00 0 2.375 0.4025 14.00
## 6 2011-04-15 00:00:00 0 2.373 0.4450 15.25
The following functions first assign event ID's to each storm event. See the gist (link above) for the assign.events function
df <- assign.events(df, value.name = "PRECIP")
print(head(df[which(!is.na(df$EVENT_ID)) - 2, ], 8)) # example showing the start of an event on row 3
## DATETIME PRECIP STAGE VELOCITY FLOW EVENT_ID EVENT_HOURS
## 48 2011-04-16 18:00:00 0.00 1.762 0.4825 13.25 NA 0
## 49 2011-04-16 19:00:00 0.00 1.782 0.4875 13.25 NA 0
## 50 2011-04-16 20:00:00 0.07 1.825 0.5325 15.00 1 1
## 51 2011-04-16 21:00:00 0.05 1.955 0.5775 17.00 1 2
## 52 2011-04-16 22:00:00 0.18 2.095 0.5475 17.00 1 3
## 53 2011-04-16 23:00:00 0.22 2.220 0.5425 17.75 1 4
## 54 2011-04-17 00:00:00 0.10 2.223 0.9550 31.50 1 5
## 55 2011-04-17 01:00:00 0.22 2.310 0.7875 27.25 1 6
## DRY_ID DRY_HOURS
## 48 1 48
## 49 1 49
## 50 NA 0
## 51 NA 0
## 52 NA 0
## 53 NA 0
## 54 NA 0
## 55 NA 0
events <- events.summary(df, value.name = "PRECIP")
head(events)
## EVENT_ID DURATION START END TOTAL PEAK
## 1 1 9 2011-04-16 20:00:00 2011-04-17 04:00:00 1.03 0.22
## 2 2 13 2011-04-19 07:00:00 2011-04-19 19:00:00 0.15 0.04
## 3 3 11 2011-04-23 02:00:00 2011-04-23 12:00:00 0.40 0.10
## 4 4 15 2011-05-07 01:00:00 2011-05-07 15:00:00 0.36 0.12
## 5 5 32 2011-05-14 17:00:00 2011-05-16 00:00:00 1.33 0.28
## 6 6 34 2011-05-16 13:00:00 2011-05-17 22:00:00 0.32 0.08
This function takes a row from the events dataframe and plots the corresponding hourly precipitation and streamflow data. The precipitation plot shows the hourly rainfall intensities and also the cumulative precipitation. The start_offset argument defines how many hours before the start of the event to start the plot. And the end_offset argument defines how many hours after the end of the event the plot should show.
plot.event <- function(event, start_offset = 0, end_offset = 48) {
x <- subset(df, DATETIME >= (event$START - hours(start_offset)) & DATETIME <=
(event$END + hours(end_offset)))
x$CUMPRECIP <- cumsum(x$PRECIP)
if (sum(is.na(x$FLOW)) < nrow(x)) {
g.pcp <- ggplot(x, aes(DATETIME, PRECIP)) + geom_bar(stat = "identity") +
geom_line(aes(y = CUMPRECIP)) + geom_vline(xint = as.numeric(event$START),
color = "red") + geom_vline(xint = as.numeric(event$END), color = "red") +
labs(x = "", y = "Precip (in/hr)")
g.q <- ggplot(x, aes(DATETIME, FLOW)) + geom_line() + geom_point() +
labs(x = "", y = "Flow (cfs)") + geom_vline(xint = as.numeric(event$START),
color = "red") + geom_vline(xint = as.numeric(event$END), color = "red")
grid.arrange(g.pcp, g.q, main = paste0("ID: ", event$EVENT_ID, ", START: ",
event$START, ", END: ", event$END, ", \nDURATION: ", event$DURATION,
" hr, PEAK: ", event$PEAK, " in/hr, TOTAL: ", event$TOTAL, " in"))
}
}
Here is an example plot for the 30th event:
events[30, ]
## EVENT_ID DURATION START END TOTAL PEAK
## 30 30 3 2011-08-21 22:00:00 2011-08-22 0.59 0.38
plot.event(events[30, ], start_offset = 6, end_offset = 48)
Note that the vertical red lines indicate the start and end datetimes of the event.
Finally, generate a PDF containing plots for all the events with total precipitation exceeding 1 in and that do not start in Jan or Feb (to avoid snow effects).
events <- subset(events, TOTAL >= 1 & !(month(START) %in% c(1, 2)))
pdf("events.pdf", width = 11, height = 8.5)
for (i in 1:nrow(event.1in)) {
plot.event(events[i, ], start_offset = 6, end_offset = 72)
}
dev.off()