Multiple log-normal probability plots in 1 plot using R
<p>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 <a href="https://stackoverflow.com/questions/11237632/log-probability-plot-in-r">here</a> to create the curve. Here is my code and a screenshot of the curve:</p>

<pre><code>library(waterData)
library(lattice)
library(latticeExtra)

site.no &lt;- "16068000" # streamflow station number
site.file &lt;- importDVs(site.no, code = "00060", stat = "00003") # reads data from NWIS
PrbGrd &lt;- qnorm(c(.001, .01, .10, .30, .50, .70, .90, .99, .999)) # create probablity scale lines
PrbGrdL&lt;-c("99.9", "99", "90", "70", "50", "30", "10", "1", "0.1") # exceedence probability scale labels
ValGrd&lt;-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&lt;-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)
       }
)
</code></pre>

<p><a href="https://i.stack.imgur.com/i9kKv.png" rel="nofollow noreferrer">Click to see my flow-duration curve</a> </p>

<p>I am having trouble plotting additional flow-duration curves on this same plot. I tried <code>lines()</code> and <code>points()</code> 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!</p>

<p>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:</p>

<p>Stations
<br>16071500
<br>16097500
<br>16017000
<br>16068000</p>

<p>The part of the code with <code>importDVs()</code> is that part that I do not know how to incorporate into <code>qqmath()</code> so that I can plot all 4 curves in one plot. Hack-R's suggestion to use <code>qqmath(~ val + val2 ...)</code> will work if I have a few number of curves to plot. For example, I could do the following:</p>

<pre><code>input.file &lt;- "Input_FDC.txt" # text file with list of station numbers
path.R &lt;- "C:/Users/ccheng/Documents/R/" # enter directory of text file
stn.file &lt;- read.delim(paste(path.R, input.file, sep = ""), header = TRUE)
site.no &lt;- 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 &lt;- importDVs(site.no[i], code = "00060", stat = "00003")
  assign(paste("flow",i,sep=""), x.1$val)
}

PrbGrd &lt;- qnorm(c(.001, .01, .10, .30, .50, .70, .90, .99, .999)) # create probablity scale lines
PrbGrdL&lt;-c("99.9", "99", "90", "70", "50", "30", "10", "1", "0.1") # exceedence probability scale labels
ValGrd&lt;-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&lt;-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
)
</code></pre>

<p>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 <code>~ flow1 + flow2 + ... flow30</code> for 30+ stations in the code. Any input is much appreciated. Thank you!  </p>