Long time no see.
Well, to feed this blog more often I decided to document some unpublished work I did in the past and also tiny things what could make others' life easier.
Around 3 years ago I wanted to plot hydrographs in R. In my case, I wanted a way that I could automatically generate a sub-hourly hydrograph to show in a very small time scale the rainfall, Q and Qbase.
I was inspired by the graphs in this very interesting paper: Aich, V., Zimmermann, A., Elsenbeer, H., 2014. Quantification and interpretation of suspended-sediment discharge hysteresis patterns: How much data do we need? Catena 122, 120–129. doi:10.1016/j.catena.2014.06.020
I found a package in R what did exactly what I wanted. The name of the package is EcoHydRology. This package is very nice, but the hydrograph could not be generated in a time scale smaller than 1 day, as I remember. To adapt it to a sub-hourly time step I just did few (and not optimal, I was still at my first steps in R) changes in the function responsible to generate those graphs. Here we go:
I call my dataset
rgraph <- 1="" 2="" function="" input="matrix(ncol" nrow="2)," precip="NULL, </p" streamflow2="NULL," streamflow="input[," timeseries="input[,"> begin = 1, endindex = length(streamflow), P.units = "mm", S.units = "",
S1.col = "black", S2.col = "red", stream.label = "Q (l/s)",
streamflow3=NULL,streamflow4=NULL, precip2=NULL)
{
if (is.null(streamflow2) & (ncol(input) > 3))
streamflow2 <- 4="" input="" p=""> if (is.null(precip) & (ncol(input) > 2)) {
precip <- 2="" input="" p=""> streamflow <- 3="" input="" p=""> }
if (!is.null(precip)) {
par(mar = c(3, 5, 1, 4))
#Here I took out some examples there was before in the lines 37-40
plot(precip, yaxt = NULL, ylab = stream.label, cex.lab = 1.2, lwd = 2,
ylim = rev(c(0, 4 * max(na.omit(precip[begin:endindex])))),
xaxt = "n", type = "S", xaxt = "n", axes = FALSE)
axis(side = 3, pos = 0, tck = 0,xaxt = "n")
axis(side = 4, lwd = 2, cex.axis = 1.3, at = seq(0, floor(max(na.omit(precip[begin:endindex])) +
1), length = (1 + ifelse(floor(max(na.omit(precip[begin:endindex])) +
1) < 10, floor(max(na.omit(precip[begin:endindex])) + 1),
4))), labels = as.integer(seq(0, floor(max(na.omit(precip[begin:endindex])) +
1), length = (1 + ifelse(floor(max(na.omit(precip[begin:endindex])) +
1) < 10, floor(max(na.omit(precip[begin:endindex])) + 1), 4)))))
if (P.units=="") {
mtext(paste("Rainfall", P.units), 4, line = 1, cex = 1.2, adj = 1)
} else mtext(paste("Rainfall (", P.units, ")", sep=""), 4, line = 2.2, cex = 1.2, adj = 1)
par(new = T)
}
if (!is.null(precip2)){
barplot(precip2[begin:endindex], yaxt = "n", space = NULL, col="pink",
ylim = rev(c(0, 4 * max(na.omit(precip[begin:endindex])))), lwd =3,
xaxt = "n")
par(new = T)
}
plot(streamflow[begin:endindex], col = S1.col, type = "l",
lwd = 3, ylab = stream.label, cex.lab = 1.2, xaxt = "n", xlab = "Time",
ylim = c(0, 1.4 * max(na.omit(streamflow[begin:endindex]), na.omit(streamflow2[begin:endindex]))),
axes = FALSE)
#mtext (expression(paste(" ", " (" , m^3/s, ")", sep="")), 2,3)
if (S.units=="m3/s" | S.units=="m3s"){
mtext (expression(paste(" (" , m^3/s, ")", sep="")), 2,1.5)
} else if (S.units=="ft3/s" | S.units=="ft3s") {
mtext (expression(paste(" (" , ft^3/s, ")", sep="")), 2,1.5)
} else if (S.units!="") mtext (paste(" (" , S.units, ")", sep=""), 2,1.5)
lines(streamflow2[begin:endindex], col = S2.col, lwd = 2,
lty = 2, xaxt = "n")
if (!is.null(streamflow3)){
lines(streamflow3[begin:endindex], col = "blue", lwd = 1, ##potential for more streamflows
lty = 3, xaxt = "n")
}
if (!is.null(streamflow4)){
lines(streamflow4[begin:endindex], col = "green", lwd = 1, ##potential for more streamflows
lty = 4, xaxt = "n")
}
axis(side = 1, lwd.ticks = 2, at = seq(1, (endindex - begin + 1), length = 14),
pos = 0, lwd = 1, cex.axis = 1.3, labels = format(timeSeries[seq(begin, endindex,
length = 14)], "%H:%M"))
axis(side = 2, pos = 0, lwd = 1, lwd.ticks = 2, cex.axis = 1.3)
legend("right", legend, c("Total Q","Baseflow"), cex = 1.2, lty=c(1,2), lwd=c(3,2), col=c("black","red"))
}
Then load my dataset and run the function:
<- 1="" 2="" function="" input="matrix(ncol" nrow="2)," precip="NULL, </p" streamflow2="NULL," streamflow="input[," timeseries="input[,"><- 4="" input="" p=""><- 2="" input="" p=""><- 3="" input="" p="">
ce1 arrow Ev2Ce (arrow bc the html makes some problems with the symbols)
ce1$time arrow strptime(ce1$time, "%H:%M")
rgraph(ce1, P.units = "mm/10 min", streamflow = ce1$Qmm, streamflow2 = ce1$Qbmm, precip = ce1$Paveg, stream.label = "Q (mm/10 min)")
text(locator(1), "P (mm) = 27.3\nImax (mm/10 min)= 11.5\nQp (mm/10 min)= 0.06\nLag time = 20 min\nRR = 0.74%")
<- ev2ce="" p=""><- ce1="" p="" strptime="" time=""> <- ev2ce="" p=""><- ce1="" p="" strptime="" time=""> <- 1="" 2="" function="" input="matrix(ncol" nrow="2)," precip="NULL, </p" streamflow2="NULL," streamflow="input[," timeseries="input[,"><- 4="" input="" p=""><- 2="" input="" p=""><- 3="" input="" p=""><- ev2ce="" p=""><- ce1="" p="" strptime="" time="">
I finally I've got this graph:
Well, to feed this blog more often I decided to document some unpublished work I did in the past and also tiny things what could make others' life easier.
Around 3 years ago I wanted to plot hydrographs in R. In my case, I wanted a way that I could automatically generate a sub-hourly hydrograph to show in a very small time scale the rainfall, Q and Qbase.
I was inspired by the graphs in this very interesting paper: Aich, V., Zimmermann, A., Elsenbeer, H., 2014. Quantification and interpretation of suspended-sediment discharge hysteresis patterns: How much data do we need? Catena 122, 120–129. doi:10.1016/j.catena.2014.06.020
I found a package in R what did exactly what I wanted. The name of the package is EcoHydRology. This package is very nice, but the hydrograph could not be generated in a time scale smaller than 1 day, as I remember. To adapt it to a sub-hourly time step I just did few (and not optimal, I was still at my first steps in R) changes in the function responsible to generate those graphs. Here we go:
I call my dataset
rgraph <- 1="" 2="" function="" input="matrix(ncol" nrow="2)," precip="NULL, </p" streamflow2="NULL," streamflow="input[," timeseries="input[,"> begin = 1, endindex = length(streamflow), P.units = "mm", S.units = "",
S1.col = "black", S2.col = "red", stream.label = "Q (l/s)",
streamflow3=NULL,streamflow4=NULL, precip2=NULL)
{
if (is.null(streamflow2) & (ncol(input) > 3))
streamflow2 <- 4="" input="" p=""> if (is.null(precip) & (ncol(input) > 2)) {
precip <- 2="" input="" p=""> streamflow <- 3="" input="" p=""> }
if (!is.null(precip)) {
par(mar = c(3, 5, 1, 4))
#Here I took out some examples there was before in the lines 37-40
plot(precip, yaxt = NULL, ylab = stream.label, cex.lab = 1.2, lwd = 2,
ylim = rev(c(0, 4 * max(na.omit(precip[begin:endindex])))),
xaxt = "n", type = "S", xaxt = "n", axes = FALSE)
axis(side = 3, pos = 0, tck = 0,xaxt = "n")
axis(side = 4, lwd = 2, cex.axis = 1.3, at = seq(0, floor(max(na.omit(precip[begin:endindex])) +
1), length = (1 + ifelse(floor(max(na.omit(precip[begin:endindex])) +
1) < 10, floor(max(na.omit(precip[begin:endindex])) + 1),
4))), labels = as.integer(seq(0, floor(max(na.omit(precip[begin:endindex])) +
1), length = (1 + ifelse(floor(max(na.omit(precip[begin:endindex])) +
1) < 10, floor(max(na.omit(precip[begin:endindex])) + 1), 4)))))
if (P.units=="") {
mtext(paste("Rainfall", P.units), 4, line = 1, cex = 1.2, adj = 1)
} else mtext(paste("Rainfall (", P.units, ")", sep=""), 4, line = 2.2, cex = 1.2, adj = 1)
par(new = T)
}
if (!is.null(precip2)){
barplot(precip2[begin:endindex], yaxt = "n", space = NULL, col="pink",
ylim = rev(c(0, 4 * max(na.omit(precip[begin:endindex])))), lwd =3,
xaxt = "n")
par(new = T)
}
plot(streamflow[begin:endindex], col = S1.col, type = "l",
lwd = 3, ylab = stream.label, cex.lab = 1.2, xaxt = "n", xlab = "Time",
ylim = c(0, 1.4 * max(na.omit(streamflow[begin:endindex]), na.omit(streamflow2[begin:endindex]))),
axes = FALSE)
#mtext (expression(paste(" ", " (" , m^3/s, ")", sep="")), 2,3)
if (S.units=="m3/s" | S.units=="m3s"){
mtext (expression(paste(" (" , m^3/s, ")", sep="")), 2,1.5)
} else if (S.units=="ft3/s" | S.units=="ft3s") {
mtext (expression(paste(" (" , ft^3/s, ")", sep="")), 2,1.5)
} else if (S.units!="") mtext (paste(" (" , S.units, ")", sep=""), 2,1.5)
lines(streamflow2[begin:endindex], col = S2.col, lwd = 2,
lty = 2, xaxt = "n")
if (!is.null(streamflow3)){
lines(streamflow3[begin:endindex], col = "blue", lwd = 1, ##potential for more streamflows
lty = 3, xaxt = "n")
}
if (!is.null(streamflow4)){
lines(streamflow4[begin:endindex], col = "green", lwd = 1, ##potential for more streamflows
lty = 4, xaxt = "n")
}
axis(side = 1, lwd.ticks = 2, at = seq(1, (endindex - begin + 1), length = 14),
pos = 0, lwd = 1, cex.axis = 1.3, labels = format(timeSeries[seq(begin, endindex,
length = 14)], "%H:%M"))
axis(side = 2, pos = 0, lwd = 1, lwd.ticks = 2, cex.axis = 1.3)
legend("right", legend, c("Total Q","Baseflow"), cex = 1.2, lty=c(1,2), lwd=c(3,2), col=c("black","red"))
}
Then load my dataset and run the function:
<- 1="" 2="" function="" input="matrix(ncol" nrow="2)," precip="NULL, </p" streamflow2="NULL," streamflow="input[," timeseries="input[,"><- 4="" input="" p=""><- 2="" input="" p=""><- 3="" input="" p="">
ce1 arrow Ev2Ce (arrow bc the html makes some problems with the symbols)
ce1$time arrow strptime(ce1$time, "%H:%M")
rgraph(ce1, P.units = "mm/10 min", streamflow = ce1$Qmm, streamflow2 = ce1$Qbmm, precip = ce1$Paveg, stream.label = "Q (mm/10 min)")
text(locator(1), "P (mm) = 27.3\nImax (mm/10 min)= 11.5\nQp (mm/10 min)= 0.06\nLag time = 20 min\nRR = 0.74%")
<- ev2ce="" p=""><- ce1="" p="" strptime="" time=""> <- ev2ce="" p=""><- ce1="" p="" strptime="" time=""> <- 1="" 2="" function="" input="matrix(ncol" nrow="2)," precip="NULL, </p" streamflow2="NULL," streamflow="input[," timeseries="input[,"><- 4="" input="" p=""><- 2="" input="" p=""><- 3="" input="" p=""><- ev2ce="" p=""><- ce1="" p="" strptime="" time="">
I finally I've got this graph:
Note that I do not use the baseflow filter of the package. I prefer to use the two parameter baseflow filter by Eckhard. For such a time scale, I do not think that any available automatic baseflow filter that I know will run that efficient.
