我的时间序列代码中存在数据空白问题。我正在从 USGS 站点提取数据,我想使用季节性平均值绘制 1986-2019 年的数据。然而,在这个特定的站点,2006-2008 年的数据差距阻止了这一点。该代码将绘制 2008-2019 年,但是当它看到数据差距时,它无法处理它。我不知道这是否像某处的 NA=TRUE 一样简单或更复杂。我只想忽略这个差距。谢谢你的帮助!
凯尔西
library(waterData)
## List of start dates
start.dates <- c ( "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",
"2017-12-01", "2018-12-01")
## List of end dates
end.dates <- c( "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",
"2018-05-31", "2019-05-31")
## Bind start and end dates into a single dataframe
flow.dates <- as.data.frame(cbind(start.dates, end.dates))
## Empty dataframe for discharge data; nrow=no of days in each inverval,
#ncol=no invtervals
discharge.data = as.data.frame(data.frame(matrix(0, nrow = 182, ncol = 33)))
## Append data for all dates
for(x in 1:33){
bin.1 <- importDVs(staid = "02196000", code = "00060", stat = "00003",
sdate = flow.dates[x, 1],
edate = flow.dates[x, 2])
bin.2 <- cleanUp(bin.1, task = "fix", replace = 0.001)
bin.3 <- fillMiss(bin.1, block = 2, pmiss = 5, model = "trend",
smooth = FALSE, log = "y")
bin.4 <- bin.3[, c(2, 3)]
colnames(bin.4) <- c(flow.dates[x, 3])
discharge.data[, x] <- bin.4[, 1]
}
## Calculate season mean & convert from ft3 to m3
mean.discharge.ws <- sapply(discharge.data, mean)
met.shift.ws <- (mean.discharge.ws / 35.315)
## Create dataframe
met.winterspring.Year <- c( 1987,
1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996,
1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005,
2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014,
2015, 2016, 2017, 2018, 2019)
df=data.frame(met.shift.ws, met.winterspring.Year)
my.formula = df$met.shift.ws ~ df$met.winterspring.Year
fit = lm(met.winterspring.Year ~ met.shift.ws ,data=df )
summary(fit)
library(ggplot2)
library(ggpmisc) # for dealing with stat equations
library (scales)
ggplot(df, aes(met.winterspring.Year, met.shift.ws)) +
geom_line(group=1) +
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() +
scale_x_continuous(breaks = breaks_width(2)) +
labs(x="", y=(expression(Discharge~(m^{3}~s^{-1}))))+
annotate(geom = "text", x = 2017, y =75, label = "Winter-Spring", hjust = "left")
dput(df[sample(1:nrow(df), 33),])
structure(list(met.shift.ws = c(23.5738634860821, 19.1108749667436,
9.29831671938426, 7.14075279970926, 4.01155689843216, 23.0833176451186,
23.3716022049591, 4.97132712961681, 7.01824863512532, 7.67068129378762,
11.5399396016697, 1.00388808416559, 4.77445533370778, 28.4664549665258,
32.4367972392891, 10.595815058508, 15.7866485772475, 6.95898919146831,
16.6857777646394, 6.97617828865174, 23.2273743529584, 5.15243499244632,
14.9389715480612, 15.3223189100295, 17.9382418515931, 4.07121541915539,
17.895891451038, 6.76626219596629, 10.9585473283619, 4.78861362338638,
31.9130961067815, 20.7926152850406, 9.99793071150851), met.winterspring.Year = c(1987,
2019, 2005, 2017, 2004, 2016, 1991, 2008, 2000, 2001, 1994, 2012,
2011, 2003, 1998, 2015, 1989, 2013, 1990, 2018, 1995, 2006, 2009,
2014, 1997, 2002, 2010, 1999, 2007, 1988, 1993, 1996, 1992)), row.names = c("X1",
"X33", "X19", "X31", "X18", "X30", "X5", "X22", "X14", "X15",
"X8", "X26", "X25", "X17", "X12", "X29", "X3", "X27", "X4", "X32",
"X9", "X20", "X23", "X28", "X11", "X16", "X24", "X13", "X21",
"X2", "X7", "X10", "X6"), class = "data.frame")
head(bin.1)
staid val dates qualcode
1 02196000 299 2018-12-01 A
2 02196000 4730 2018-12-02 A
3 02196000 1810 2018-12-03 A
4 02196000 554 2018-12-04 A
5 02196000 283 2018-12-05 A
6 02196000 187 2018-12-06 A
>
> head(discharge.data)
X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 X16 X17 X18 X19
1 1480 14 181 93 114 33 395 82 711 145 529 1000 25 20 40 4.83 16.3 34 98.0
2 1780 22 136 88 103 36 329 62 374 135 1860 459 26 16 33 3.87 15.8 41 86.3
3 1220 25 111 85 101 43 277 53 276 125 694 270 26 13 30 3.00 14.7 40 90.8
4 583 19 97 82 112 168 238 48 274 120 298 206 26 12 26 2.53 14.6 73 81.2
5 342 17 88 80 144 239 223 134 2050 122 218 187 27 12 24 2.22 17.4 122 69.1
6 244 15 79 80 130 147 212 293 873 120 184 168 27 12 24 2.18 21.2 95 64.3
X20 X21 X22 X23 X24 X25 X26 X27 X28 X29 X30 X31 X32 X33
1 17.8 41.0 6.58 1120.0 112 26.7 31.1 0.43 35.3 27.4 134 53.6 24.1 299
2 17.5 45.4 6.82 360.0 738 68.1 26.1 0.41 29.6 21.5 138 55.6 28.6 4730
3 18.0 44.1 9.69 149.0 5210 59.3 18.1 0.40 27.9 16.8 242 54.3 31.7 1810
4 18.5 35.6 8.90 83.8 1540 33.0 13.5 0.47 50.1 13.2 217 33.1 32.7 554
5 25.2 31.1 7.47 58.2 467 23.9 11.1 0.56 74.1 11.2 156 27.6 33.2 283
6 680.0 28.4 6.73 45.0 292 19.0 12.5 0.52 79.0 11.3 129 27.1 35.8 187