Create 3 new vignettes

Using this as a template: https://code.usgs.gov/water/wrtdsplus/-/blob/main/vignettes/Introduction.Rmd

Create 3 vignettes that we'll use to describe the 3 different workflows. For now, we won't have anything substantial in them. As we find appropriate data sets, we'll make them more complete.

For the vignette that uses the "WRTZS" (2nd workflow) we could start working on setting up the same data that Matt uses here:

https://code.usgs.gov/mdiebel/egret-hs

So here's how you get the data:

library(EGRET)
site = "05427930"
pcode = "00665"
start_date = "2012-10-01"
end_date = "2021-09-30"
INFO = readNWISInfo(site, pcode, interactive=FALSE)
INFO$shortName = "Dorn Creek"
INFO$staAbbrev = "DCM"
INFO$paramShortName = "Total Phosphorus"
INFO$constitAbbrev = "TP"
INFO$windowY = 7
INFO$windowS = 0.25
INFO$windowQ = 2
INFO$minNumObs = 100

Sample <- readNWISSample(site, pcode, start_date, end_date)
Daily <- readNWISDaily(site, "00060", start_date, end_date)

eList <- mergeReport(INFO, Daily, Sample)
plot(eList)

Then, to calculate "Z", you can use the BaseflowSeparation function:

# BaseflowSeparation function from EcoHydRology package

BaseflowSeparation <-
  function(streamflow, filter_parameter=0.925, passes=3){
    suppressWarnings(Ends<-c(1,length(streamflow))*rep(1,(passes+1))) # Start and end values for the filter function
    suppressWarnings(AddToStart<-c(1,-1)*rep(1,passes))
    btP<-streamflow##Previous pass's baseflow approximation
    qft<-vector(length=length(streamflow))
    bt<-vector(length=length(streamflow))
    bt[1]<-if(streamflow[1]<quantile(streamflow,0.25)) streamflow[1] else mean(streamflow)/1.5
    ##Guess baseflow value in first time step.  
    for(j in 1:passes){
      for (i in (Ends[j]+AddToStart[j]):Ends[j+1]){
        if ((filter_parameter*bt[i-AddToStart[j]]+((1-filter_parameter)/2)*(btP[i]+btP[i-AddToStart[j]]))>btP[i]){
          bt[i]<-btP[i]
        } else bt[i]<-filter_parameter*bt[i-AddToStart[j]]+((1-filter_parameter)/2)*(btP[i]+btP[i-AddToStart[j]])
        qft[i]<-streamflow[i]-bt[i]
      }
      if (j<passes){
        btP<-bt
        bt[Ends[j+1]]<-if(streamflow[Ends[j+1]]<mean(btP))streamflow[Ends[j+1]]/1.2 else mean(btP)
        ##Refines the approximation of end values after the first pass
      }
    }
    f <- data.frame(bt,qft)
    return(f)
  }

And here's the script to use the log of quickflow for Z:

Qs <- BaseflowSeparation(Daily$Q, 
                         filter_parameter = 0.925,
                         passes = 3)

# to make an offset for when quickflow is 0:
Qq05 <- quantile(Qs$qft[Qs$qft > 0], probs = 0.05)

Daily$Qb <- Qs$bt
Daily$LogQb <- log(Daily$Qb)
Daily$Qq <- Qs$qft
Daily$LogQq <- log(Qs$qft + Qq05)

Daily$Z <- Daily$LogQq

Sample <- dplyr::left_join(Sample, 
                           Daily[, c("Date", "Z", "Qb", "Qq")],
                           by = "Date")

eList <- mergeReport(INFO, Daily, Sample)
plot(eList)

library(WRTDSplus)
eList <- modelEstimation_ReplaceQ(eList)
plotConcHist(eList)