Данные временных рядов. Можно ли агрегировать сезоны? Хочу заговор только зима-весна

Question

Мне интересно, может ли кто-нибудь помочь мне найти способ агрегирования сезонов. Я хочу построить только зимне-весенний сезон для разгрузки данных за 1920-2019 годы. Я пытался использовать dplyr, но он позволяет только для года, месяца или среднесуточного значения. Любая помощь будет принята с благодарностью.

Ниже приведен мой код для среднего значения дня.

library(dataRetrieval)
siteNo <- "02202500"
pCode <- "00060"

daily <- readNWISdv(siteNo, pCode, "1920-10-01","2019-09-30")
daily <- renameNWISColumns(daily)
head(daily)

attach(daily)
library(dplyr)
library(ggpmisc) 
library(ggplot2)

df= data.frame(Date, Flow)
Flow2 = (0.0283168)*Flow
my.formula <- Flow2 ~ Date


p1 <- ggplot(data = df,
             aes(Date, Flow2)) +
  geom_line()+
  geom_smooth(method = "lm", se=FALSE) +
  stat_poly_eq(formula = my.formula,
               eq.with.lhs = "italic(hat(y))~`=`~",
               aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               parse = TRUE)+
  theme_classic()+
  labs(x="", y=(expression(Discharge~(m^{3}~s^{-1})))) 
p1

Ниже приведена моя попытка кода найти среднее время года

library(dataRetrieval)
library(plyr)
library(dplyr)
library(lubridate) # for working with dates
siteNo <- "02202500"
pCode <- "00060"
start.dates <- c("1969-12-01", "1970-12-01", "1971-12-01", "1972-12-01",
                 "1973-12-01", "1974-12-01", "1975-12-01", "1976-12-01",
                 "1977-12-01", "1978-12-01", "1979-12-01", "1980-12-01",
                 "1981-12-01", "1982-12-01", "1983-12-01", "1984-12-01",
                 "1985-12-01", "1986-12-01", "1987-12-01", "1988-12-01",
                 "1989-12-01", "1990-12-01", "1991-12-01", "1992-12-01",
                 "1993-12-01", "1994-12-01", "1995-12-01", "1996-12-01",
                 "1997-12-01", "1998-12-01", "1999-12-01", "2000-12-01",
                 "2001-12-01", "2002-12-01", "2003-12-01", "2004-12-01",
                 "2005-12-01", "2006-12-01", "2007-12-01", "2008-12-01",
                 "2009-12-01", "2010-12-01", "2011-12-01", "2012-12-01",
                 "2013-12-01", "2014-12-01", "2015-12-01", "2016-12-01")


## List of end dates
end.dates <- c("1970-05-31", "1971-05-31", "1972-05-30", "1973-05-31", 
               "1974-05-31", "1975-05-31", "1976-05-30", "1977-05-31", 
               "1978-05-31", "1979-05-31", "1980-05-30", "1981-05-31", 
               "1982-05-31", "1983-05-31", "1984-05-30", "1985-05-31", 
               "1986-05-31", "1987-05-31", "1988-05-30", "1989-05-31", 
               "1990-05-31", "1991-05-31", "1992-05-30", "1993-05-31", 
               "1994-05-31", "1995-05-31", "1996-05-30", "1997-05-31", 
               "1998-05-31", "1999-05-31", "2000-05-30", "2001-05-31", 
               "2002-05-31", "2003-05-31", "2004-05-30", "2005-05-31", 
               "2006-05-31", "2007-05-31", "2008-05-30", "2009-05-31", 
               "2010-05-31", "2011-05-31", "2012-05-30", "2013-05-31", 
               "2014-05-31", "2015-05-31", "2016-05-30", "2017-05-31")
daily = readNWISdv(siteNo, pCode, start.dates, end.dates, statCd="00003")
daily = renameNWISColumns(daily)
dates= format(as.Date(daily$Date), format="%Y-%m-%d")
dates = as.data.frame(cbind(start.dates, end.dates))
mean_Flow = sapply(daily, mean)
mean_shift = (mean_Flow/ 35.315)

library(ggpmisc) # for dealing with stat equations
library(ggplot2) # for making plots 
library(scales) #for working with date_format

df$dates = as.Date(df$dates, format = "%Y-%m-%d")
my.formula = df$mean_shift ~ df$dates

p1=ggplot(df,aes(dates, mean_shift)) +
  geom_line() +
  geom_smooth(method = "lm", se=FALSE) +
  stat_poly_eq(formula = my.formula,
               eq.with.lhs = "italic(hat(y))~`=`~",
               aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), parse = TRUE)+
  theme_classic()+
  labs(x="", y=(expression(Discharge~(m^{3}~s^{-1}))))+
  scale_x_date(breaks = "5 year", labels = date_format("%Y-%m"))+
  scale_y_continuous(limits = c(0,600)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
p1