Isolating peaks in from time series data in R

Question

I am working with 15-min streamflow data (one column of datetime, one of streamflow). I would like to write code that finds peaks (storms) over a certain threshold and select the data after until the data drops back down to a baseflow threshold. I would like to isolate these peak-and-recession sections, and then select those same time periods from other time-series data of water quality metrics.

I have tried the findpeaks function of the pracma package, and the output it gives me is the index of the peak and its start and endpoints, but I don't know how to convert these indices back to datetimes so I can select from the other time series.

I also have tried the find.spell.lengths from the Hystrostats package. It returns the spell length but I don't know how to convert that back to datetime intervals. Also this package seems to be designed for daily data (not 15-minute data) so I don't know if that will mess up the results.

Is there a different package or function I should be using, or a better way to do this?

Answer 1

I think you are almost there. I just simulated some data and use the pracma package. You mentioned having problems getting the dates back, what you do is use the index returned from the findpeak call, and subset your dataframe according to that,

For example, to get datetime of the absolute peak, do:

peaks <-findpeaks(stream_data$streamflow, minpeakheight = 1, minpeakdistance = 4, sortstr = FALSE)
stream_data[peaks[,2],"datetime"]

so see below code on how to call them back.

library(pracma)
#simulate data using pracma findpeaks example
x <- seq(0, 1, len = 1024)
     pos <- c(0.1, 0.13, 0.15, 0.23, 0.25, 0.40, 0.44, 0.65, 0.76, 0.78, 0.81)
     hgt <- c(4, 5, 3, 4, 5, 4.2, 2.1, 4.3, 3.1, 5.1, 4.2)
     wdt <- c(0.005, 0.005, 0.006, 0.01, 0.01, 0.03, 0.01, 0.01, 0.005, 0.008, 0.005)

pSignal <- numeric(length(x))
for (i in seq(along=pos)) {
             pSignal <- pSignal + hgt[i]/(1 + abs((x - pos[i])/wdt[i]))^4
     }

#give a datetime
START = as.POSIXlt("2018-09-07 10:00:00")
END = as.POSIXlt("2018-09-16 10:00:00")
#make your data frame
stream_data = data.frame(
       datetime=seq(START,END,length.out=length(pSignal)),
       streamflow=pSignal)
#plot
with(stream_data,plot(datetime,streamflow,type="l"))
#call peaks like you did
peak <-findpeaks(stream_data$streamflow, minpeakheight = 1, minpeakdistance = 4, sortstr = FALSE)
## call
data.frame(
       peak_number = 1:nrow(peak),
       peak_max = stream_data[peak[,2],"datetime"],
       peak_heigth = stream_data[peak[,2],"streamflow"],
       peak_start = stream_data[peak[,3],"datetime"],
       peak_end = stream_data[peak[,4],"datetime"]
)