Newer
Older
library(zoo)
library(StatsAnalysisTrend)
# Sourcing R file
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
# Compute the time gap by station
get_lacune = function (df_data, df_info) {
# Get all different stations code
Code = levels(factor(df_info$code))
# Create new vector to stock results for cumulative time gap by station
tLac = c()
# Create new vector to stock results for mean time gap by station
meanLac = c()
# Get rows where there is no NA
NoNA = complete.cases(df_data)
# Get data where there is no NA
df_data_NoNA = df_data[NoNA,]
# For every station
for (code in Code) {
# Get only the data rows for the selected station
df_data_code = df_data[df_data$code==code,]
# Get date for the selected station
Date = df_data_code$Date
# Get time span for the selection station
span = as.numeric(Date[length(Date)] - Date[1])
# Get only the data rows with no NA for the selected station
df_data_NoNA_code = df_data_NoNA[df_data_NoNA$code==code,]
# Get date for the selected station
Date_NoNA = df_data_NoNA_code$Date
# Compute the time gap
lac = as.numeric(diff(Date_NoNA) - 1)
# Compute the cumulative gap
lac_sum = sum(lac)
# Store the cumulative gap rate
tLac = c(tLac, lac_sum/span)
# Compute the mean gap
lac_mean = mean(lac[lac != 0])
# Store the mean gap
meanLac = c(meanLac, lac_mean)
}
# Compute the cumulative gap rate in pourcent
tLac100 = tLac * 100
# Create a tibble
df_lac = tibble(code=Code, tLac100=tLac100, meanLac=meanLac)
return (df_lac)
}
get_intercept = function (df_Xtrend, df_Xlist, unit2day=365.25) {
for (g in df_Xlist$info$group) {
df_data_code = df_Xlist$data[df_Xlist$data$group == g,]
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
104
105
106
107
108
109
110
111
112
113
114
115
116
df_Xtrend_code = df_Xtrend[df_Xtrend$group == g,]
Start = df_Xtrend_code$period_start
UStart = levels(factor(Start))
End = df_Xtrend_code$period_end
UEnd = levels(factor(End))
nPeriod = max(length(UStart), length(UEnd))
for (i in 1:nPeriod) {
df_data_code_per =
df_data_code[df_data_code$Date >= Start[i]
& df_data_code$Date <= End[i],]
df_Xtrend_code_per =
df_Xtrend_code[df_Xtrend_code$period_start == Start[i]
& df_Xtrend_code$period_end == End[i],]
id = which(df_Xtrend$group == g
& df_Xtrend$period_start == Start[i]
& df_Xtrend$period_end == End[i])
mu_X = mean(df_data_code_per$Qm3s, na.rm=TRUE)
mu_t = as.numeric(mean(c(Start[i],
End[i]),
na.rm=TRUE)) / unit2day
b = mu_X - mu_t * df_Xtrend_code_per$trend
df_Xtrend$intercept[id] = b
}
}
return (df_Xtrend)
}
get_period = function (per, df_Xtrend, df_XEx, df_Xlist) {
df_Xtrend = tibble(df_Xtrend)
df_Xtrend$period_start = as.Date("1970-01-01")
df_Xtrend$period_end = as.Date("1970-01-01")
df_Xlisttmp = reprepare(df_XEx, df_Xlist, colnamegroup=c('code'))
df_XExtmp = df_Xlisttmp$data
for (g in df_Xlisttmp$info$group) {
df_XExtmp_code = df_XExtmp[df_XExtmp$group == g,]
iStart = which.min(abs(df_XExtmp_code$Date
- as.Date(per[1])))
iEnd = which.min(abs(df_XExtmp_code$Date
- as.Date(per[2])))
id = which(df_Xtrend$group1 == g)
df_Xtrend$period_start[id] =
as.Date(df_XExtmp_code$Date[iStart])
df_Xtrend$period_end[id] =
as.Date(df_XExtmp_code$Date[iEnd])
}
get_QAtrend = function (df_data, period) {
# AVERAGE ANNUAL FLOW : QA #
period = as.list(period)
Imax = 0
df_QAtrendB = tibble()
df_QAlist = prepare(df_data, colnamegroup=c('code'))
df_QAEx = extract.Var(data.station=df_QAlist,
funct=mean,
timestep='year',
period=per,
pos.datetime=1,
na.rm=TRUE)
df_QAtrend = Estimate.stats(data.extract=df_QAEx)
I = interval(per[1], per[2])
if (I > Imax) {
Imax = I
df_QAlistB = df_QAlist
df_QAExB = df_QAEx
}
df_QAtrend = get_period(per, df_QAtrend, df_QAEx, df_QAlist)
df_QAtrendB = bind_rows(df_QAtrendB, df_QAtrend)
res_QAtrend = clean(df_QAtrendB, df_QAExB, df_QAlistB)
get_QMNAtrend = function (df_data, period) {
# MONTHLY MINIMUM FLOW IN THE YEAR : QMNA #
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
period = as.list(period)
Imax = 0
df_QMNAtrendB = tibble()
for (per in period) {
df_QMNAlist = prepare(df_data, colnamegroup=c('code'))
df_QMNAEx = extract.Var(data.station=df_QMNAlist,
funct=mean,
period=per,
timestep='month',
pos.datetime=1,
na.rm=TRUE)
df_QMNAlist = reprepare(df_QMNAEx, df_QMNAlist, colnamegroup=c('code'))
df_QMNAEx = extract.Var(data.station=df_QMNAlist,
funct=min,
period=per,
timestep='year',
pos.datetime=1,
na.rm=TRUE)
df_QMNAtrend = Estimate.stats(data.extract=df_QMNAEx)
I = interval(per[1], per[2])
if (I > Imax) {
Imax = I
df_QMNAlistB = df_QMNAlist
df_QMNAExB = df_QMNAEx
}
df_QMNAtrend = get_period(per, df_QMNAtrend, df_QMNAEx,
df_QMNAlist)
df_QMNAtrendB = bind_rows(df_QMNAtrendB, df_QMNAtrend)
}
res_QMNAtrend = clean(df_QMNAtrendB, df_QMNAExB, df_QMNAlistB)
get_VCN10trend = function (df_data, df_meta, period) {
# MINIMUM 10 DAY AVERAGE FLOW OVER THE YEAR : VCN10 #
# Get all different stations code
df_data_roll = tibble()
for (c in Code) {
df_data_code = df_data[df_data$code == c,]
df_data_code = tibble(Date=rollmean(df_data_code$Date,
10,
fill=NA),
Qm3s=rollmean(df_data_code$Qm3s,
10,
fill=NA),
code=c)
df_data_roll = bind_rows(df_data_roll, df_data_code)
}
period = as.list(period)
Imax = 0
df_VCN10trendB = tibble()
for (per in period) {
df_VCN10list = prepare(df_data_roll, colnamegroup=c('code'))
df_VCN10Ex = extract.Var(data.station=df_VCN10list,
funct=min,
period=per,
timestep='year',
pos.datetime=1,
na.rm=TRUE)
df_VCN10trend = Estimate.stats(data.extract=df_VCN10Ex)
I = interval(per[1], per[2])
if (I > Imax) {
Imax = I
df_VCN10listB = df_VCN10list
df_VCN10ExB = df_VCN10Ex
}
df_VCN10trend = get_period(per, df_VCN10trend, df_VCN10Ex,
df_VCN10list)
df_VCN10trendB = bind_rows(df_VCN10trendB, df_VCN10trend)
}
res_VCN10trend = clean(df_VCN10trendB, df_VCN10ExB, df_VCN10listB)