# Usefull library
library(dplyr)
library(zoo)
library(StatsAnalysisTrend)
library(lubridate)
# Sourcing R file
source('processing/format.R', encoding='latin1')
# 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) {
df_Xtrend$intercept = NA
for (g in df_Xlist$info$group) {
df_data_code = df_Xlist$data[df_Xlist$data$group == g,]
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])
}
return (df_Xtrend)
}
# get_QAtrend = function (df_data, period, p_thresold) {
# # AVERAGE ANNUAL FLOW : QA #
# period = as.list(period)
# for (p_thr in p_thresold) {
# Imax = 0
# df_QAtrendB = tibble()
# for (per in period) {
# 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)
# return (res_QAtrend)
# }
get_QAtrend = function (df_data, period, p_thresold) {
# AVERAGE ANNUAL FLOW : QA #
period = as.list(period)
Imax = 0
df_QAtrendB = tibble()
for (per in period) {
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,
level=p_thresold)
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)
return (res_QAtrend)
}
get_QMNAtrend = function (df_data, period, p_thresold) {
# MONTHLY MINIMUM FLOW IN THE YEAR : QMNA #
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_QMNAEx = extract.Var(data.station=df_QMNAlist,
funct=mean,
period=per,
timestep='year-month',
per.start="01",
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,
level=p_thresold)
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)
return (res_QMNAtrend)
}
get_VCN10trend = function (df_data, df_meta, period, p_thresold) {
# MINIMUM 10 DAY AVERAGE FLOW OVER THE YEAR : VCN10 #
# Get all different stations code
Code = levels(factor(df_meta$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,
level=p_thresold)
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)
return (res_VCN10trend)
}