r - How to Create a Graph of Statistical Time Series

Question

I have data in the following format:

        Date    Year    Month   Day     Flow
1   1953-10-01  1953    10       1      530
2   1953-10-02  1953    10       2      530
3   1953-10-03  1953    10       3      530

I would like to create a graph like this:

Here is my current image and code:

library(ggplot2)
library(plyr)
library(reshape2)
library(scales)

## Read Data
df <- read.csv("Salt River Flow.csv")

## Convert Date column to R-recognized dates
df$Date <- as.Date(df$Date, "%m/%d/%Y")

## Finds Water Years (Oct - Sept)
df$WY <- as.POSIXlt(as.POSIXlt(df$Date)+7948800)$year+1900

## Normalizes Water Years so stats can be applied to just months and days
df$w <- ifelse(month(df$Date) %in% c(10,11,12), 1903, 1904)

##Creates New Date (dat) Column
df$dat <- as.Date(paste(df$w,month(df$Date),day(df$Date), sep = "-"))

## Creates new data frame with summarised data by MonthDay
PlotData <- ddply(df, .(dat), summarise, Min = min(Flow), Tenth = quantile(Flow, p = 0.05), TwentyFifth = quantile(Flow, p =    0.25), Median = quantile(Flow, p = 0.50), Mean = mean(Flow), SeventyFifth = quantile(Flow, p = 0.75), Ninetieth = quantile(Flow, p = 0.90), Max = max(Flow))

## Melts data so it can be plotted with ggplot
m <- melt(PlotData, id="dat")

## Plots
p <- ggplot(m, aes(x = dat)) + 
geom_ribbon(aes(min = TwentyFifth, max = Median), data = PlotData, fill = alpha("black", 0.1), color = NA) + 
geom_ribbon(aes(min = Median, max = SeventyFifth), data = PlotData, fill = alpha("black", 0.5), color = NA) + 
scale_x_date(labels = date_format("%b"), breaks = date_breaks("month"), expand = c(0,0)) + 
geom_line(data = subset(m, variable == "Mean"), aes(y = value), size = 1.2) + 
theme_bw() + 
geom_line(data = subset(m, variable %in% c("Min","Max")), aes(y = value, group = variable)) + 
geom_line(data = subset(m, variable %in% c("Ninetieth","Tenth")), aes(y = value, group = variable), linetype = 2) + 
labs(x = "Water Year", y = "Flow (cfs)")

p

I am very close but there are some issues I'm having. First, if you can see a way to improve my code, please let me know. The main problem I ran into was that I needed two dataframes to make this graph: one melted, and one not. The unmelted dataframe was necessary (I think) to create the ribbons. I tried many ways to use the melted dataframe for the ribbons, but there was always a problem with the aesthetic length.

Second, I know to have a legend - and I want one, I need to have something in the aesthetics of each line/ribbon, but I am having trouble getting that to work. I think it would involve scale_fill_manual.

Third, and I don't know if this is possible, I would like to have each month label in between the tick marks, not on them (like in the above image).

Any help is greatly appreciated (especially with creating more efficient code).

Thank you.

Answer 1

Something along these lines might get you close with base:

library(lubridate)
library(reshape2)
# simulating data...
Date  <- seq(as.Date("1953-10-01"),as.Date("2010-10-01"),by="day")
Year  <- year(Date)
Month <- month(Date)
Day <- day(Date)
set.seed(1)
Flow <- rpois(length(Date), 2000)
Data <- data.frame(Date=Date,Year=Year,Month=Month,Day=Day,Flow=Flow)

# use acast to get it in a convenient shape:
PlotData <- acast(Data,Year~Month+Day,value.var="Flow")
# apply for quantiles
Quantiles <- apply(PlotData,2,function(x){
    quantile(x,probs=c(1,.9,.75,.5,.25,.1,0),na.rm=TRUE)
  })
Mean <- colMeans(PlotData, na.rm=TRUE)
# ugly way to get month tick separators
MonthTicks <- cumsum(table(unlist(lapply(strsplit(names(Mean),split="_"),"[[",1))))

# and finally your question:
plot(1:366,seq(0,max(Flow),length=366),type="n",xlab = "Water Year",ylab="Discharge",axes=FALSE)
polygon(c(1:366,366:1),c(Quantiles["50%",],rev(Quantiles["75%",])),border=NA,col=gray(.6))
polygon(c(1:366,366:1),c(Quantiles["50%",],rev(Quantiles["25%",])),border=NA,col=gray(.4))
lines(1:366,Quantiles["90%",], col = gray(.5), lty=4)
lines(1:366,Quantiles["10%",], col = gray(.5))
lines(1:366,Quantiles["100%",], col = gray(.7))
lines(1:366,Quantiles["0%",], col = gray(.7), lty=4)
lines(1:366,Mean,lwd=3)
axis(1,at=MonthTicks, labels=NA)
text(MonthTicks-15,-100,1:12,pos=1,xpd=TRUE)
axis(2)

The plotting code really isn't that tricky. You'll need to clean up the aesthetics, but polygon() is usually my strategy for shaded regions in plots (confidence bands, whatever).