-
Notifications
You must be signed in to change notification settings - Fork 0
/
locate_prec_event.Rmd
104 lines (86 loc) · 3.03 KB
/
locate_prec_event.Rmd
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
96
97
98
99
100
101
102
103
```{r}
library('tidyverse')
library('glue')
library(plotly)
## Global variables ##
PLACE <- 'Miedes'
VWC_THRESHOLD <- 0.03
THRESHOLD_CHANGE <- 0.0075
COMPARISON_WINDOW <- 4 * 4 # four hours
OUTPUT_PATH <- 'output-data'
```
```{r}
db.env <- read.csv("output-data/clim-allsites.csv")
db.env <- db.env %>% select(ts, site, vwc) %>% mutate(ts = as_datetime(ts, tz = 'Europe/Madrid'))
head(db.env)
```
```{r}
db.env <- db.env %>% filter(site == PLACE)
```
```{r}
# 1) locate intervals where there has been a significant rain episode (after 4 hours, there was a significant increase of VWC)
prec_changes_bools <- diff(db.env$vwc, lag=COMPARISON_WINDOW) >= VWC_THRESHOLD
# prec_changes <- db.env[prec_changes_bools,]
# 2) Set each interval of a prec event (consecutive timestamps) in a different item of a list of events
# prec_events <- list()
prec_events_pos <- list()
i = 1;
l_idx = 1;
while( i <= length(prec_changes_bools)) {
if (prec_changes_bools[i]) {
# prec_events[[l_idx]] <- db.env[i,]
prec_events_pos[[l_idx]] <- i
l_idx = l_idx + 1;
while(prec_changes_bools[i+1]) { i = i + 1 }
}
i = i + 1;
}
# while ( consecutive(inst, inst + 1) ) {
# prev_events [[i]].append( inst ++ )
# }
```
```{r}
# 3) locale the changing point in each precipitation event (first change from current trend to new trend of precipitation)
diff.next <- diff(db.env$vwc, lag=1)
i = 1
changing_point.df <- data.frame()
for (pos in prec_events_pos) {
interval <- diff.next[pos:(pos+COMPARISON_WINDOW)] # interval of comparison window
change <- interval > THRESHOLD_CHANGE
changing_point_pos <- pos + first (which ( change ) ) - 1 # we have to substract one to get the moment where the change starts
if(is.na(changing_point_pos)) {
r <- data.frame(ts = db.env[pos,]$ts, site = db.env[pos,]$site, vwc = NA)
changing_point.df <- rbind(changing_point.df,r)
} else {
changing_point.df <- rbind(changing_point.df,db.env[changing_point_pos,])
}
i = i + 1;
}
```
```{r}
print(changing_point.df)
```
^ fill the rest of the missing prec events start (NAs) by hand:
```{r}
missing.prec_events <- c() %>%
map(\(dt) as.POSIXct(dt, tz = 'Europe/Madrid'))
missing.prec_events.df <- subset(db.env, ts %in% missing.prec_events)
changing_point.df <- rbind.data.frame(changing_point.df, missing.prec_events.df)
```
Clean NAs (missing starting prec events, introduced by hand), duplicated (long raining episodes), and sort by ts.
```{r}
changing_point.df <- drop_na(changing_point.df, vwc)
changing_point.df <- changing_point.df[!duplicated(changing_point.df), ]
changing_point.df <- changing_point.df %>% arrange(ts)
print(changing_point.df)
```
Visually, double check
```{r}
fig <- plot_ly() %>%
add_trace(data = db.env, x = ~ts, y = ~vwc, type = 'scatter', mode = 'lines', name = 'TMS') %>%
add_trace(data = changing_point.df, x = ~ts, y = ~vwc, type = 'scatter', mode = 'markers', name = 'Prec events')
fig
```
```{r}
write_csv(changing_point.df, file.path(OUTPUT_PATH, glue('prec_events-{PLACE}.csv')), append = F, col_names = T)
```