diff --git a/processing/analyse.R b/processing/analyse.R
index 851f63e9ee2a90be636850a5dc2ea9814d262ad4..5977ded39318d64477ed4fb5a3799142611132f6 100644
--- a/processing/analyse.R
+++ b/processing/analyse.R
@@ -1,5 +1,10 @@
# Usefull library
library(dplyr)
+library(zoo)
+library(StatsAnalysisTrend)
+
+# Sourcing R file
+source('processing/format.R')
# Compute the time gap by station
@@ -57,3 +62,139 @@ get_lacune = function (df_data, df_info) {
return (df_lac)
}
+
+
+get_intercept = function (df_Xtrend, df_Xlist, unit2day=365.25) {
+
+ intercept = c()
+ # Group = levels(factor())
+ for (g in df_Xlist$info$group) {
+ df_data_code = df_Xlist$data[df_Xlist$data$group == g,]
+
+ trend = df_Xtrend$trend[df_Xtrend$group == g]
+ mu_X = mean(df_data_code$Qm3s, na.rm=TRUE)
+
+ mu_t = as.numeric(mean(df_data_code$Date,
+ na.rm=TRUE))/unit2day
+
+ b = mu_X - mu_t * trend
+
+ intercept = append(intercept, b)
+ }
+ return (intercept)
+}
+
+
+get_QAtrend = function (df_data, period) {
+ # AVERAGE ANNUAL FLOW : QA #
+ ### /!\ verify order conservation ###
+ df_QAlist = prepare(df_data, colnamegroup=c('code'))
+
+ df_QAEx = extract.Var(data.station=df_QAlist,
+ funct=mean,
+ timestep='year',
+ period=period,
+ pos.datetime=1,
+ na.rm=TRUE)
+
+ df_QAtrend = Estimate.stats(data.extract=df_QAEx)
+
+ res_QAtrend = clean(df_QAtrend, df_QAEx, df_QAlist)
+
+ return (res_QAtrend)
+}
+
+get_QMNAtrend = function (df_data, period) {
+ # MONTHLY MINIMUM FLOW IN THE YEAR : QMNA #
+ df_QMNAlist = prepare(df_data, colnamegroup=c('code'))
+
+ ### /!\ PLUS RAPIDE ###
+ # fMNA = function (X) {
+ # # prendre un paquet de 1 ans et faire la moyenne par mois et retourner le minimum des debit
+ # dpm = length(X)/12
+ # # print(dpm)
+ # # print(length(X))
+ # monthmean = c()
+ # for (i in 1:12) {
+ # id = round(dpm*(i-1)+1, 0)
+ # iu = round(i*dpm, 0)
+ # monthmean = append(monthmean, mean(X[id:iu], na.rm=TRUE))
+ # # print(paste('start', id))
+ # # print(paste('end', iu))
+ # # print('')
+ # }
+ # # print(monthmean)
+ # return (min(monthmean, na.rm=TRUE))
+ # }
+ # df_QMNAEx = extract.Var(data.station=df_QMNAlist,
+ # funct=fMNA,
+ # period=period,
+ # pos.datetime=1)#,
+ # na.rm=TRUE) ### /!\ PAS COMPRIS ###
+
+ df_QMNAEx = extract.Var(data.station=df_QMNAlist,
+ funct=mean,
+ period=period,
+ timestep='month',
+ pos.datetime=1,
+ na.rm=TRUE)
+
+ ### /!\ NOM DE COLONNE PAS CONSERVER ###
+ df_QMNAlist = reprepare(df_QMNAEx, df_QMNAlist, colnamegroup=c('code'))
+
+ df_QMNAEx = extract.Var(data.station=df_QMNAlist,
+ funct=min,
+ period=period,
+ timestep='year',
+ pos.datetime=1,
+ na.rm=TRUE)
+
+ df_QMNAtrend = Estimate.stats(data.extract=df_QMNAEx)
+
+ res_QMNAtrend = clean(df_QMNAtrend, df_QMNAEx, df_QMNAlist)
+
+ return (res_QAtrend)
+}
+
+
+get_VCN10trend = function (df_data, df_info, period) {
+ # MINIMUM 10 DAY AVERAGE FLOW OVER THE YEAR : VCN10 #
+
+ # Get all different stations code
+ Code = levels(factor(df_info$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)
+
+ }
+
+ df_VCN10list = prepare(df_data_roll, colnamegroup=c('code'))
+
+ df_VCN10Ex = extract.Var(data.station=df_VCN10list,
+ funct=min,
+ period=period,
+ timestep='year',
+ pos.datetime=1,
+ na.rm=TRUE)
+
+ df_VCN10trend = Estimate.stats(data.extract=df_VCN10Ex)
+
+ res_VCN10trend = clean(df_VCN10trend, df_VCN10Ex, df_VCN10list)
+
+ return (res_VCN10trend)
+}
+
+
+
diff --git a/processing/format.R b/processing/format.R
index e050c26dc2634676d6e09b50b4b5b58cbea828cb..28f9c7ee1559188d301f60da7292aab255012e40 100644
--- a/processing/format.R
+++ b/processing/format.R
@@ -50,41 +50,49 @@ prepare = function(df_data, colnamegroup=NULL) {
}
-reprepare = function(df_dataEx, df_list, colnamegroup=NULL) {
+reprepare = function(df_XEx, df_Xlist, colnamegroup=NULL) {
- colnames(df_dataEx) = c('Date', 'group', 'Qm3s')
+ colnames(df_XEx) = c('Date', 'group', 'Qm3s')
- df_dataEx$Date = as.character(df_dataEx$Date)
- exDate = df_dataEx$Date[1]
+ df_XEx$Date = as.character(df_XEx$Date)
+ exDate = df_XEx$Date[1]
nbt = lengths(regmatches(exDate, gregexpr('-', exDate)))
if (nbt == 1) {
- df_dataEx$Date = paste(df_dataEx$Date, '01', sep='-')
+ df_XEx$Date = paste(df_XEx$Date, '01', sep='-')
} else if (nbt == 0) {
- df_dataEx$Date = paste(df_dataEx$Date, '01', '01', sep='-')
+ df_XEx$Date = paste(df_XEx$Date, '01', '01', sep='-')
} else if (nbt != 2) {
stop('erreur of date format')
}
- df_dataEx = bind_cols(Date=as.Date(df_dataEx$Date,
+ df_XEx = bind_cols(Date=as.Date(df_XEx$Date,
format="%Y-%m-%d"),
- df_dataEx[-1],
- df_list$info[df_dataEx$group,
- 2:ncol(df_list$info)])
+ df_XEx[-1],
+ df_Xlist$info[df_XEx$group,
+ 2:ncol(df_Xlist$info)])
- df_listEx = prepare(df_dataEx, colnamegroup=colnamegroup)
- return (df_listEx)
+ df_XlistEx = prepare(df_XEx, colnamegroup=colnamegroup)
+ return (df_XlistEx)
}
+clean = function (df_Xtrend, df_XEx, df_Xlist) {
-clean = function (df_xTrend, df_list) {
+ # print(str(df_XEx))
- df_xTrend = bind_cols(df_xTrend,
- df_list$info[df_xTrend$group1,
- 2:ncol(df_list$info)])
+ df_Xlist = reprepare(df_XEx, df_Xlist, colnamegroup=c('code'))
- colnames(df_xTrend)[1] = 'group'
+ df_Xtrend = bind_cols(df_Xtrend,
+ df_Xlist$info[df_Xtrend$group1,
+ 2:ncol(df_Xlist$info)])
- return (df_xTrend)
+ colnames(df_Xtrend)[1] = 'group'
+
+ intercept = get_intercept(df_Xtrend, df_Xlist, unit2day=365.25)
+
+ df_Xtrend$intercept = intercept
+ df_Xtrend = relocate(df_Xtrend, intercept, .after=trend)
+
+ return (list(trend=df_Xtrend, data=df_Xlist$data, info=df_Xlist$info))
}
diff --git a/script.R b/script.R
index 723373c222cb142eaf7f7d4185ad9e20cd986ea0..1ff1e56a7e90409833f7f1f4286ef0d267b26ea0 100644
--- a/script.R
+++ b/script.R
@@ -79,8 +79,7 @@ source('processing/analyse.R')
source('plotting/panel.R')
# Usefull library
-library(StatsAnalysisTrend)
-library(dplyr)
+
# Result directory
resdir = file.path(computer_work_path, 'results')
@@ -150,63 +149,11 @@ df_info = df_join$info
# df_lac = get_lacune(df_data, df_info)
-# AVERAGE ANNUAL FLOW : QA #
-### /!\ verify order conservation ###
-df_QAlist = prepare(df_data, colnamegroup=c('code'))
-
-df_QAEx = extract.Var(data.station=df_QAlist,
- funct=mean,
- timestep='year',
- period=period,
- pos.datetime=1,
- na.rm=TRUE)
-df_QATrend = Estimate.stats(data.extract=df_QAEx)
-df_QATrend = clean(df_QATrend, df_QAlist)
-
-
-# MONTHLY MINIMUM FLOW IN THE YEAR : QMNA #
-df_QMNAlist = prepare(df_data, colnamegroup=c('code'))
-
-### /!\ PLUS RAPIDE ###
-# fMNA = function (X) {
-# # prendre un paquet de 1 ans et faire la moyenne par mois et retourner le minimum des debit
-# dpm = length(X)/12
-# # print(dpm)
-# # print(length(X))
-# monthmean = c()
-# for (i in 1:12) {
-# id = round(dpm*(i-1)+1, 0)
-# iu = round(i*dpm, 0)
-# monthmean = append(monthmean, mean(X[id:iu], na.rm=TRUE))
-# # print(paste('start', id))
-# # print(paste('end', iu))
-# # print('')
-# }
-# # print(monthmean)
-# return (min(monthmean, na.rm=TRUE))
-# }
-# df_QMNAEx = extract.Var(data.station=df_QMNAlist,
-# funct=fMNA,
-# period=period,
-# pos.datetime=1)#,
- # na.rm=TRUE) ### /!\ PAS COMPRIS ###
-
-df_QMNAEx = extract.Var(data.station=df_QMNAlist,
- funct=mean,
- period=period,
- timestep='month',
- pos.datetime=1,
- na.rm=TRUE)
-
-### /!\ NOM DE COLONNE PAS CONSERVER ###
-df_QMNAlist = reprepare(df_QMNAEx, df_QMNAlist, colnamegroup=c('code'))
-
-df_QMNAEx = extract.Var(data.station=df_QMNAlist,
- funct=min,
- period=period,
- timestep='year',
- pos.datetime=1,
- na.rm=TRUE)
-
-df_QMNATrend = Estimate.stats(data.extract=df_QMNAEx)
-df_QMNATrend = clean(df_QMNATrend, df_QMNAlist)
+# QA TREND #
+res_QAtrend = get_QAtrend(df_data, period)
+
+# QMNA TREND #
+# res_QMNAtrend = get_QMNAtrend(df_data, period)
+
+# VCN10 TREND #
+# res_VCN10trend = get_VCN10trend(df_data, df_info, period)
diff --git a/script_install.R b/script_install.R
index dbd980c206766c55074d27713b605eca3a8fbdc0..b1383666477110d139b05c656a23e944bf9071fd 100644
--- a/script_install.R
+++ b/script_install.R
@@ -6,6 +6,7 @@ install.packages("dplyr")
install.packages("ggplot2")
install.packages("officer")
install.packages("lubridate")
+install.packages('zoo')
library(devtools)
install_github("https://github.com/benRenard/BFunk") #type '1'