2

I am a newbie working with streamflow duration curves and the function fdc. I am working with more than 300 series and I am interested in saving the low quartile threshold Qlow.thr value that appears in the plot generated:

enter image description here

Here is the reproducible example:

dat <- c(13.05, 90.29,  5.68, 49.13, 26.39, 15.06, 23.39, 17.98,  4.21,  2.51, 38.29,  8.57,  2.48 , 3.78, 18.09 ,15.16, 13.46,  8.69, 6.85, 11.97, 12.10,  9.87 ,21.89,  2.60  ,2.40, 27.40,  4.94, 83.17 ,12.10,  5.08 ,12.42,  6.19  ,3.60 ,32.58, 53.69, 38.49,3.61, 14.84, 34.48,  1.91, 21.79, 31.53,  6.70,  9.52, 22.64,  1.80 , 8.13, 10.60, 12.73,  4.17,  6.70 ,16.45)

fdc(dat,plot = T,lQ.thr=0.8,ylab='Hm3',main='Upstream monthly duration curve',thr.shw=TRUE)

The fdc function returns a vector of probabilities, but I am not sure how to convert these probabilities to the original units and select the 80% percentile value expressed in Hm3 as I would do with pnorm, for example, in case of working with normal probabilities.

Thank you so much.

Marina
  • 143
  • 1
  • 10
  • Related https://stackoverflow.com/questions/59914776/flow-duration-curve-using-facet-wrap-of-ggplot-in-r – Tung Feb 04 '20 at 01:33

1 Answers1

2

You can construct the FDC yourself by 

dat <- c(13.05, 90.29,  5.68, 49.13, 26.39, 15.06, 23.39, 17.98,  
         4.21,  2.51, 38.29,  8.57,  2.48 , 3.78, 18.09 ,15.16,
         13.46,  8.69, 6.85, 11.97, 12.10,  9.87 ,21.89,  2.60,
         2.40, 27.40,  4.94, 83.17 ,12.10,  5.08 ,12.42,  6.19,
         3.60 ,32.58, 53.69, 38.49,3.61, 14.84, 34.48,  1.91, 
         21.79, 31.53,  6.70,  9.52, 22.64,  1.80 , 8.13, 10.60, 
         12.73,  4.17,  6.70 ,16.45)

dat <- sort(dat, decreasing = T)
df  <- data.frame(x = 100/length(dat) * 1:length(dat), y = dat)

plot(x = df$x, y = df$y, type = "l", log = "y")

The plot created using <code>fdc</code> versus the one using <code>plot</code>

So the sorted flow data is simply plotted against the percentage exceedance scale. This scale is created by dividing 100% by the number of data points which gives us the increment for each point. Therefore

quantile(dat, p = c(0.2, 0.8), type = 1)

gives you your desired results.

Notice that the computation of the quantile differs in fdc. It seems like they just use

p <- c(0.8, 0.2)
dat[round(p * length(dat))]

> [1]  4.21 27.40

to compute the values.

Martin Schmelzer
  • 23,283
  • 6
  • 73
  • 98
  • Way simpler than what I got to do (modifying the original function to return the lQ). Although the q20 (4.356) slightly varies from the value given by the FDC function (4.21), the solution is great. Thank you so so much Martin. – Marina Oct 16 '18 at 09:51
  • 2
    I updated my answer. There are many ways to compute quantiles. The actual way used in `fdc` is added at the bottom of my answer. – Martin Schmelzer Oct 16 '18 at 09:53
  • Very clear and didactic, thank you so much once again Martin. – Marina Oct 16 '18 at 09:58