-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathexample.R
95 lines (69 loc) · 3.52 KB
/
example.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
rm(list = ls())
"%>%" = magrittr::`%>%`
#### hobo dataset ####
hobo <- read.csv("./data/hobo_quincemil.csv", header = F) %>%
setNames(c("date_time", "count")) %>%
transform(date_time = strptime(date_time, "%m/%d/%Y %H:%M:%S")) %>%
transform(diff = c(0, diff(date_time, 1))) %>%
# Deleting (equal to 0) FALSE DOUBLES (differences in count less or equal to 1) and
# LOGGED (NA values in count) date_times
transform(pp = ifelse(diff <= 1 | is.na(count), 0, 0.2))
# Deleting (equal to 0) events in days(hours) of INSTALLATION and DATA DOWNLOADING
# This corresponds to the first and two last events, respectively
hobo[c(1, 13918, 13919), "pp"] = 0
## hobo dataset to continuous date (to seconds and minute)
# as no events is equal to 0, you could create a continuous date-time
# with all non numeric values to 0 OR interpolate the numeric values
# into the NA values (of the continuous date_time)
# from minute
hobo_minute <- xts::xts(hobo$pp, hobo$date_time) %>%
xts::period.apply(x = .,
INDEX = xts::endpoints(., on = "minutes"),
FUN = sum)
# making seconds to 00
zoo::index(hobo_minute) = as.POSIXct(paste(format(time(hobo_minute), "%Y-%m-%d %H:%M"), ":00", sep = ""))
# creating a continuous time span
xts_minute = seq(as.POSIXct("2017-02-14 00:00:00"), as.POSIXct("2017-07-24 23:00:00"), by = 60)
# converting irregular to regular xts ts
hobo_minute = merge(hobo_minute, xts_minute)
hobo_minute[is.na(hobo_minute)] = 0
# deleting (setting NA) values outside of the range
hobo_minute[ time(hobo_minute) < hobo$date_time[1]] = NA
hobo_minute[ time(hobo_minute) > hobo$date_time[length(hobo$date_time)]] = NA
## hobo dataset to continuous date (to hours and daily)
hobo_hourly <- xts::xts(hobo$pp, hobo$date_time) %>%
xts::period.apply(x = .,
INDEX = xts::endpoints(., on = "hours"),
FUN = sum)
zoo::index(hobo_hourly) = as.POSIXct(paste(format(time(hobo_hourly), "%Y-%m-%d %H"), ":00:00", sep = ""))
xts_hourly = seq(as.POSIXct("2017-02-14 00:00:00"), as.POSIXct("2017-07-24 23:00:00"), by = 3600)
hobo_hourly = merge(hobo_hourly, xts_hourly)
hobo_hourly[is.na(hobo_hourly)] = 0
hobo_hourly[ time(hobo_hourly) < hobo$date_time[1]] = NA
hobo_hourly[ time(hobo_hourly) > hobo$date_time[length(hobo$date_time)]] = NA
# to daily, it's just needed minute or hour object, in this case the first option is used
hobo_daily <- hobo_hourly %>%
xts::lag.xts(., k = -8) %>% # daily sum as in a conventional station (PERU time 7am-7pm)
xts::apply.daily(sum, na.rm = FALSE)
# converting time to date
zoo::index(hobo_daily) = as.Date(format(time(hobo_daily), "%Y-%m-%d"))
#### conventional dataset ######
conventional <- read.csv("./data/conventional_quincemil.csv", header = F, na.strings = -999) %>%
setNames(c("Y", "M", "D", "pp")) %>%
transform(date_time = ISOdate(year = Y, month = M, day = D))
conventional_daily = xts::xts(conventional$pp, as.Date(conventional$date_time))
#### comparison ######
cbind(hobo_daily, conventional_daily) %>%
window(start = "2017-02-01", end = "2017-07-31") %>%
xts::plot.xts(col = c("royalblue", "black"), lwd = c(2, 3.5), lty = c(1, 1))
xts::addLegend("top", on = 1,
legend.names = c("hobo_daily", "conventional_daily"),
lwd = c(2, 3.5),
col = c("royalblue", "black"),
cex = 1)
cor(cbind(hobo_daily, conventional_daily) %>% zoo::coredata(),
use = "pairwise.complete.obs",
method = "spearman")
(hobo_daily-conventional_daily) %>%
.[!is.na(.)] %>%
mean()