diff --git a/plotting/layout.R b/plotting/layout.R
index b2b8d361fdb671f6d1c10878a871e09983f77b2d..f2e27d0d125eeb4f5bddbee85b9816344bde5130 100644
--- a/plotting/layout.R
+++ b/plotting/layout.R
@@ -91,7 +91,6 @@ panels_layout = function (df_data, df_meta, layout_matrix, figdir='', filedir_op
dir.create(outdirTmp)
}
-
# Get all different stations code
Code = levels(factor(df_meta$code))
nCode = length(Code)
diff --git a/plotting/panel.R b/plotting/panel.R
index 9a9ce89f0a07f9c31069a815c3d56a2b5a42e5b1..d3ae1d1695c6fbad2c604d5247541ed0ecdd6a84 100644
--- a/plotting/panel.R
+++ b/plotting/panel.R
@@ -10,8 +10,7 @@ library(ggh4x)
library(RColorBrewer)
-time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missRect=FALSE, unit2day=365.25, period=NULL, last=FALSE, color=NULL, norm=FALSE) {
-
+time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missRect=FALSE, unit2day=365.25, period=NULL, last=FALSE, color=NULL) {
if (type == 'sqrt(Q)') {
df_data_code$Qm3s = sqrt(df_data_code$Qm3s)
@@ -20,43 +19,21 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR
maxQ = max(df_data_code$Qm3s, na.rm=TRUE)
power = get_power(maxQ)
- if (norm) {
- dbrk = 10^power
- } else {
- dbrk = 1
- }
-
- df_data_code$Qm3sN = df_data_code$Qm3s / dbrk
-
- if (!is.null(df_trend_code)) {
- df_trend_code$trendN = df_trend_code$trend / dbrk
- df_trend_code$interceptN = df_trend_code$intercept / dbrk
- }
-
- maxQN = max(df_data_code$Qm3sN, na.rm=TRUE)
-
maxQtmp = maxQ/10^power
if (maxQtmp >= 5) {
dbrk = 1.0
- accuracy = 0.1
} else if (maxQtmp < 5 & maxQtmp >= 3) {
dbrk = 0.5
- accuracy = 0.1
} else if (maxQtmp < 3 & maxQtmp >= 2) {
dbrk = 0.4
- accuracy = 0.1
} else if (maxQtmp < 2 & maxQtmp >= 1) {
dbrk = 0.2
- accuracy = 0.1
} else if (maxQtmp < 1) {
dbrk = 0.1
- accuracy = 0.1
}
- if (!norm) {
- dbrk = dbrk * 10^power
- accuracy = NULL
- }
+ dbrk = dbrk * 10^power
+ accuracy = NULL
dDate = as.numeric(df_data_code$Date[length(df_data_code$Date)] -
df_data_code$Date[1]) / unit2day
@@ -118,17 +95,17 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR
if (type == 'sqrt(Q)' | type == 'Q') {
p = p +
- geom_line(aes(x=df_data_code$Date, y=df_data_code$Qm3sN),
+ geom_line(aes(x=df_data_code$Date, y=df_data_code$Qm3s),
color='grey20',
size=0.3)
} else {
p = p +
- geom_point(aes(x=df_data_code$Date, y=df_data_code$Qm3sN),
+ geom_point(aes(x=df_data_code$Date, y=df_data_code$Qm3s),
shape=1, color='grey20', size=1)
}
if (missRect) {
- NAdate = df_data_code$Date[is.na(df_data_code$Qm3sN)]
+ NAdate = df_data_code$Date[is.na(df_data_code$Qm3s)]
dNAdate = diff(NAdate)
NAdate_Down = NAdate[append(Inf, dNAdate) != 1]
NAdate_Up = NAdate[append(dNAdate, Inf) != 1]
@@ -137,7 +114,7 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR
geom_rect(aes(xmin=NAdate_Down,
ymin=0,
xmax=NAdate_Up,
- ymax=maxQN*1.1),
+ ymax=maxQ*1.1),
linetype=0, fill='Wheat', alpha=0.4)
}
@@ -147,74 +124,124 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR
geom_rect(aes(xmin=min(df_data_code$Date),
ymin=0,
xmax=period[1],
- ymax= maxQN*1.1),
+ ymax= maxQ*1.1),
linetype=0, fill='grey85', alpha=0.3) +
geom_rect(aes(xmin=period[2],
ymin=0,
xmax=max(df_data_code$Date),
- ymax= maxQN*1.1),
+ ymax= maxQ*1.1),
linetype=0, fill='grey85', alpha=0.3)
}
+
if (!is.null(df_trend_code)) {
- if (df_trend_code$p <= p_threshold) {
+
+ # print(df_trend_code)
- abs = c(df_data_code$Date[1],
- df_data_code$Date[length(df_data_code$Date)])
+ Start = df_trend_code$period_start
+ UStart = levels(factor(Start))
+ End = df_trend_code$period_end
+ UEnd = levels(factor(End))
+
+ nPeriod = max(length(UStart), length(UEnd))
+
+ Periods = vector(mode='list', length=nPeriod)
+ # for (i in 1:nPeriod) {
+ # Periods[[i]] = as.Date(c(Period_start[i], Period_end[i]))
+ # }
+
+ ltype = c('solid', 'dashed', 'dotted', 'twodash')
+ for (i in 1:nPeriod) {
- abs_num = as.numeric(abs) / unit2day
+ df_trend_code_per =
+ df_trend_code[df_trend_code$period_start == Start[i]
+ & df_trend_code$period_end == End[i],]
- ord = abs_num * df_trend_code$trendN +
- df_trend_code$interceptN
+ # print(df_trend_code_per)
- if (!is.null(color)) {
- p = p +
- geom_line(aes(x=abs, y=ord),
- color=color,
- size=0.7)
- } else {
- p = p +
- geom_line(aes(x=abs, y=ord),
- color='cornflowerblue')
- }
-
- if (norm) {
- p = p +
- ggtitle(bquote(bold(.(type))~~'['*m^{3}*'.'*s^{-1}*'] x'~10^{.(as.character(power))}~~~bold('tendance')~.(format(df_trend_code$trend, scientific=TRUE, digits=3))~'['*m^{3}*'.'*s^{-1}*'.'*an^{-1}*']'))
- } else {
- p = p +
- ggtitle(bquote(bold(.(type))~~'['*m^{3}*'.'*s^{-1}*']'~~~bold('tendance')~.(format(df_trend_code$trend, scientific=TRUE, digits=3))~'['*m^{3}*'.'*s^{-1}*'.'*an^{-1}*']'))
- }
-
- } else {
- if (norm) {
- p = p +
- ggtitle(bquote(bold(.(type))~~'['*m^{3}*'.'*s^{-1}*'] x'~10^{.(as.character(power))}~~~bold('tendance')~.(format(df_trend_code$trend, scientific=TRUE, digits=3))~'['*m^{3}*'.'*s^{-1}*'.'*an^{-1}*']'))
- } else {
- p = p +
- ggtitle(bquote(bold(.(type))~~'['*m^{3}*'.'*s^{-1}*']'~~~bold('tendance')~.(format(df_trend_code$trend, scientific=TRUE, digits=3))~'['*m^{3}*'.'*s^{-1}*'.'*an^{-1}*']'))
+ if (df_trend_code_per$p <= p_threshold) {
+
+ iStart = which.min(abs(df_data_code$Date - Start[i]))
+ iEnd = which.min(abs(df_data_code$Date - End[i]))
+
+ abs = c(df_data_code$Date[iStart],
+ df_data_code$Date[iEnd])
+
+ # abs = seq(df_data_code$Date[1],
+ # df_data_code$Date[length(df_data_code$Date)],
+ # length.out=10)
+
+ # abs[abs <= df_data_code$Date[iStart]] = NA
+ # abs[abs >= df_data_code$Date[iEnd]] = NA
+
+
+ # print(abs)
+ # print(df_trend_code_per$trend)
+ # print(df_trend_code_per$intercept)
+
+
+ abs_num = as.numeric(abs) / unit2day
+
+
+ ord = abs_num * df_trend_code_per$trend +
+ df_trend_code_per$intercept
+
+ plot = tibble(abs=abs, ord=ord)
+
+ if (!is.null(color)) {
+ p = p +
+ geom_line(data=plot, aes(x=abs, y=ord),
+ color=color[i],
+ linetype=ltype[i], size=0.7)
+ } else {
+ p = p +
+ geom_line(aes(x=abs, y=ord),
+ color='cornflowerblue')
+ }
}
}
- } else {
- if (norm) {
- p = p +
- ggtitle(bquote(bold(.(type))~' ['*m^{3}*'.'*s^{-1}*'] x'~10^{.(as.character(power))}))
- } else {
- p = p +
- ggtitle(bquote(bold(.(type))~' ['*m^{3}*'.'*s^{-1}*']'))
- }
}
+
+
+
+
+ # if (norm) {
+ # p = p +
+ # ggtitle(bquote(bold(.(type))~~'['*m^{3}*'.'*s^{-1}*'] x'~10^{.(as.character(power))}~~~bold('tendance')~.(format(df_trend_code$trend, scientific=TRUE, digits=3))~'['*m^{3}*'.'*s^{-1}*'.'*an^{-1}*']'))
+ # } else {
+ # p = p +
+ # ggtitle(bquote(bold(.(type))~~'['*m^{3}*'.'*s^{-1}*']'~~~bold('tendance')~.(format(df_trend_code$trend, scientific=TRUE, digits=3))~'['*m^{3}*'.'*s^{-1}*'.'*an^{-1}*']'))
+ # }
+
+ # } else {
+ # if (norm) {
+ # p = p +
+ # ggtitle(bquote(bold(.(type))~~'['*m^{3}*'.'*s^{-1}*'] x'~10^{.(as.character(power))}~~~bold('tendance')~.(format(df_trend_code$trend, scientific=TRUE, digits=3))~'['*m^{3}*'.'*s^{-1}*'.'*an^{-1}*']'))
+ # } else {
+ # p = p +
+ # ggtitle(bquote(bold(.(type))~~'['*m^{3}*'.'*s^{-1}*']'~~~bold('tendance')~.(format(df_trend_code$trend, scientific=TRUE, digits=3))~'['*m^{3}*'.'*s^{-1}*'.'*an^{-1}*']'))
+ # }
+ # }
+ # } else {
+ # if (norm) {
+ # p = p +
+ # ggtitle(bquote(bold(.(type))~' ['*m^{3}*'.'*s^{-1}*'] x'~10^{.(as.character(power))}))
+ # } else {
+ # p = p +
+ # ggtitle(bquote(bold(.(type))~' ['*m^{3}*'.'*s^{-1}*']'))
+ # }
+ # }
- if (norm) {
- p = p +
- ylab(bquote('débit ['*m^{3}*'.'*s^{-1}*'] x'~10^{.(as.character(power))}))
- } else {
- p = p +
- ylab(expression(paste('débit [', m^{3}, '.',
- s^{-1}, ']', sep='')))
- }
+ # if (norm) {
+ # p = p +
+ # ylab(bquote('débit ['*m^{3}*'.'*s^{-1}*'] x'~10^{.(as.character(power))}))
+ # } else {
+ # p = p +
+ # ylab(expression(paste('débit [', m^{3}, '.',
+ # s^{-1}, ']', sep='')))
+ # }
p = p +
# xlab('date') +
@@ -229,8 +256,8 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR
expand=c(0, 0))
p = p +
- scale_y_continuous(breaks=seq(0, maxQN*10, dbrk),
- limits=c(0, maxQN*1.1),
+ scale_y_continuous(breaks=seq(0, maxQ*10, dbrk),
+ limits=c(0, maxQ*1.1),
expand=c(0, 0),
labels=label_number(accuracy=accuracy))
diff --git a/processing/analyse.R b/processing/analyse.R
index 83d7c1aec14eab7f63b0df660d5f579fdaecbf43..953c7ec12c7b48fa390fabb47eb14537ff84cb19 100644
--- a/processing/analyse.R
+++ b/processing/analyse.R
@@ -65,23 +65,76 @@ get_lacune = function (df_data, df_info) {
get_intercept = function (df_Xtrend, df_Xlist, unit2day=365.25) {
+
+ df_Xtrend$intercept = NA
- intercept = c()
- # Group = levels(factor())
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")
- trend = df_Xtrend$trend[df_Xtrend$group == g]
- mu_X = mean(df_data_code$Qm3s, na.rm=TRUE)
+ df_Xlisttmp = reprepare(df_XEx, df_Xlist, colnamegroup=c('code'))
+ df_XExtmp = df_Xlisttmp$data
- mu_t = as.numeric(mean(df_data_code$Date,
- na.rm=TRUE))/unit2day
+ for (g in df_Xlisttmp$info$group) {
+
+ df_XExtmp_code = df_XExtmp[df_XExtmp$group == g,]
- b = mu_X - mu_t * trend
+ iStart = which.min(abs(df_XExtmp_code$Date
+ - as.Date(per[1])))
+ iEnd = which.min(abs(df_XExtmp_code$Date
+ - as.Date(per[2])))
- intercept = append(intercept, b)
+ 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 (intercept)
+ return (df_Xtrend)
}
@@ -93,7 +146,7 @@ get_QAtrend = function (df_data, period) {
Imax = 0
df_QAtrendB = tibble()
- for (per in period){
+ for (per in period) {
df_QAlist = prepare(df_data, colnamegroup=c('code'))
@@ -113,11 +166,10 @@ get_QAtrend = function (df_data, period) {
df_QAExB = df_QAEx
}
- df_QAtrend = bind_cols(df_QAtrend,
- tibble(period_start=as.Date(per[1])),
- tibble(period_end=as.Date(per[2])))
- df_QAtrendB = bind_rows(df_QAtrendB, df_QAtrend)
+ 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)
@@ -125,6 +177,7 @@ get_QAtrend = function (df_data, period) {
return (res_QAtrend)
}
+
get_QMNAtrend = function (df_data, period) {
# MONTHLY MINIMUM FLOW IN THE YEAR : QMNA #
@@ -162,9 +215,9 @@ get_QMNAtrend = function (df_data, period) {
df_QMNAExB = df_QMNAEx
}
- df_QMNAtrend = bind_cols(df_QMNAtrend,
- tibble(period_start=as.Date(per[1])),
- tibble(period_end=as.Date(per[2])))
+ df_QMNAtrend = get_period(per, df_QMNAtrend, df_QMNAEx,
+ df_QMNAlist)
+
df_QMNAtrendB = bind_rows(df_QMNAtrendB, df_QMNAtrend)
}
@@ -222,9 +275,9 @@ get_VCN10trend = function (df_data, df_meta, period) {
df_VCN10ExB = df_VCN10Ex
}
- df_VCN10trend = bind_cols(df_VCN10trend,
- tibble(period_start=as.Date(per[1])),
- tibble(period_end=as.Date(per[2])))
+ df_VCN10trend = get_period(per, df_VCN10trend, df_VCN10Ex,
+ df_VCN10list)
+
df_VCN10trendB = bind_rows(df_VCN10trendB, df_VCN10trend)
}
diff --git a/processing/format.R b/processing/format.R
index 7ac9101f8a02e2a06568e414e20a3252d0bc2f04..e60138773eedff650f86051c4dbec10f71bc0b6c 100644
--- a/processing/format.R
+++ b/processing/format.R
@@ -100,9 +100,9 @@ clean = function (df_Xtrend, df_XEx, df_Xlist) {
colnames(df_Xtrend)[1] = 'group'
- intercept = get_intercept(df_Xtrend, df_Xlist, unit2day=365.25)
+ df_Xtrend = get_intercept(df_Xtrend, df_Xlist, unit2day=365.25)
- df_Xtrend$intercept = intercept
+ # 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 c34e2d92ddfcc3828ace843cc8d4984b68d14672..d45154af561e0892c4544d8fc98b026db7029e5f 100644
--- a/script.R
+++ b/script.R
@@ -16,19 +16,20 @@ computer_work_path =
### BANQUE HYDRO ###
# Path to the directory where BH data is stored
BHfiledir =
- # ""
- "BanqueHydro_Export2021"
+ ""
+ # "BanqueHydro_Export2021"
## Manual selection ##
# Name of the file that will be analysed from the BH directory
BHfilename =
- # ""
- c("S2235610_HYDRO_QJM.txt",
- "P1712910_HYDRO_QJM.txt",
- "P0885010_HYDRO_QJM.txt",
- "A1000030_HYDRO_QJM.txt",
- "A2250310_HYDRO_QJM.txt"
- )
+ ""
+ # c(
+ # "S2235610_HYDRO_QJM.txt",
+ # "P1712910_HYDRO_QJM.txt",
+ # "P0885010_HYDRO_QJM.txt",
+ # "A1000030_HYDRO_QJM.txt",
+ # "A2250310_HYDRO_QJM.txt"
+ # )
## Or list selection ##
# Path to the list file of BH data that will be analysed
@@ -43,13 +44,13 @@ BHlistname =
### NIVALE ###
# Path to the directory where NV data is stored
NVfiledir =
- ""
- # "France207"
+ # ""
+ "France207"
# Name of the file that will be analysed from the NV directory
NVfilename =
- ""
- # "all"
+ # ""
+ "all"
# Path to the list file of metadata about station that will be analysed
@@ -57,14 +58,14 @@ NVlistdir =
""
NVlistname =
- ""
- # "liste_bv_principaux_global.txt"
+ # ""
+ "liste_bv_principaux_global.txt"
### TREND ANALYSIS ###
# Time period to analyse
-period_all = c("1700-01-01", "2020-12-31")
-period2 = c("1968-01-01", "2020-12-31")
+period_all = c("1800-01-01", "2019-12-31")
+period2 = c("1980-01-01", "2019-12-31")
########################
@@ -148,13 +149,14 @@ df_meta = df_join$meta
# QA TREND #
-# res_QAtrend = get_QAtrend(df_data, period=list(period_all, period2))
+res_QAtrend = get_QAtrend(df_data, period=list(period_all, period2))
# QMNA TREND #
-res_QMNAtrend = get_QMNAtrend(df_data, period=list(period_all, period2))
+# res_QMNAtrend = get_QMNAtrend(df_data, period=list(period_all, period2))
# VCN10 TREND #
-# res_VCN10trend = get_VCN10trend(df_data, period=list(period_all, period2))
+res_VCN10trend = get_VCN10trend(df_data, df_meta,
+ period=list(period_all, period2))
# TIME PANEL #
@@ -203,6 +205,22 @@ panels_layout(list(res_QAtrend$data, res_VCN10trend$data),
filename_opt='')
+
+# panels_layout(list(res_QAtrend$data),
+# layout_matrix=c(1),
+# df_meta=df_meta,
+# df_trend=list(res_QAtrend$trend),
+# type=list(bquote(Q[A])),
+# missRect=list(TRUE),
+# period=period_all,
+# info_header=TRUE,
+# time_header=df_data,
+# time_ratio=2,
+# var_ratio=5,
+# figdir=figdir,
+# filename_opt='')
+
+
### /!\ Removed 185 row(s) containing missing values (geom_path) -> remove NA ###