Aloha. I am trying to plot multiple flow-duration curves in one plot. A flow-duration curve shows the probability of streamflow and has log-normal y-axis and probability scale x-axis. I used the example here to create the curve. Here is my code and a screenshot of the curve:
library(waterData)
library(lattice)
library(latticeExtra)
site.no <- "16068000" # streamflow station number
site.file <- importDVs(site.no, code = "00060", stat = "00003") # reads data from NWIS
PrbGrd <- qnorm(c(.001, .01, .10, .30, .50, .70, .90, .99, .999)) # create probablity scale lines
PrbGrdL<-c("99.9", "99", "90", "70", "50", "30", "10", "1", "0.1") # exceedence probability scale labels
ValGrd<-c(seq(0.001,0.01,0.001),seq(0.01,0.1,0.01),seq(0.1,1,0.1),seq(1,10,1),seq(10,100,10)) # create value grid lines
ValGrd<-log10(ValGrd) # convert value grid lines to log
qqmath(~ val, # val = data column with streamflow data
data = site.file,
distribution = function(p) qnorm(p),
main = "Flow-duration curves for streamgages on Kauai",
pch = 20,
cex = 0.5,
xlab = "Percentage of time discharge was equaled or exceeded (%)",
ylab = "Daily mean discharge, in cubic feet per second",
xlim = c(max(PrbGrd),min(PrbGrd)),
scales = list(y=list(log=10,alternating=1),x = list(at = PrbGrd, labels = PrbGrdL, cex = 0.8)),
yscale.components=yscale.components.log10ticks,
panel = function(x,...){
panel.abline(v=PrbGrd ,col="grey",lty=3)
panel.abline(h=ValGrd,col="grey",lty=3)
panel.qqmath(x,distribution=qnorm)
}
)
Click to see my flow-duration curve
I am having trouble plotting additional flow-duration curves on this same plot. I tried lines() and points() as suggested in other posts but they do not work. Ideally, I would like to have a text file containing a list of streamflow station numbers and one plot with flow-duration curves for all the streamflow stations listed in the text file. Any input is much appreciated. Thank you!
In response to Hack-R's comment, here is an example. I have a text file named "Input_FDC.txt" with a list of station numbers. For simplicity, I have 4 station numbers (but in actuality, I have 30+ stations). Here are the contents of the text file:
Stations
16071500
16097500
16017000
16068000
The part of the code with importDVs() is that part that I do not know how to incorporate into qqmath() so that I can plot all 4 curves in one plot. Hack-R's suggestion to use qqmath(~ val + val2 ...) will work if I have a few number of curves to plot. For example, I could do the following:
input.file <- "Input_FDC.txt" # text file with list of station numbers
path.R <- "C:/Users/ccheng/Documents/R/" # enter directory of text file
stn.file <- read.delim(paste(path.R, input.file, sep = ""), header = TRUE)
site.no <- as.character(stn.file$Stations) # for importDVs to work, station no. needs to be in quotes
for(i in 1:length(site.no)) {
x.1 <- importDVs(site.no[i], code = "00060", stat = "00003")
assign(paste("flow",i,sep=""), x.1$val)
}
PrbGrd <- qnorm(c(.001, .01, .10, .30, .50, .70, .90, .99, .999)) # create probablity scale lines
PrbGrdL<-c("99.9", "99", "90", "70", "50", "30", "10", "1", "0.1") # exceedence probability scale labels
ValGrd<-c(seq(0.001,0.01,0.001),seq(0.01,0.1,0.01),seq(0.1,1,0.1),seq(1,10,1),seq(10,100,10)) # create value grid lines
ValGrd<-log10(ValGrd) # convert value grid lines to log
qqmath(~ flow1 + flow2 + flow3 + flow4,
distribution = function(p) qnorm(p),
main = "Flow-duration curves for streamgages on Kauai",
auto.key = list(points = F, lines = T, text = site.no, columns = 4),
pch = 20,
cex = 0.5,
xlab = "Percentage of time discharge was equaled or exceeded (%)",
ylab = "Daily mean discharge, in cubic feet per second",
xlim = c(max(PrbGrd),min(PrbGrd)),
scales = list(y=list(log=10,alternating=1),x = list(at = PrbGrd, labels = PrbGrdL, cex = 0.8)),
yscale.components=yscale.components.log10ticks
)
But with 30+ curves, I am looking for a solution that would allow me to loop through the list of stations in the "Input_FDC.txt" file and not having to include ~ flow1 + flow2 + ... flow30 for 30+ stations in the code. Any input is much appreciated. Thank you!
lines()andpoints()(or evenadd=T) won't work is thatqqmath()plots using thelatticepackage which does not allow new layers like this. – vincentmajor