diff --git a/plotting/layout.R b/plotting/layout.R
index f5072a77df3b5230ccd0ea1182ac989dbb40face..54d399da1d74b535bbf7790ad42c96a79c68f626 100644
--- a/plotting/layout.R
+++ b/plotting/layout.R
@@ -14,14 +14,34 @@ library(RColorBrewer)
source('plotting/panel.R', encoding='latin1')
-panels_layout = function (df_data, df_meta, layout_matrix, figdir='', filedir_opt='', filename_opt='', variable='', df_trend=NULL, p_threshold=0.1, unit2day=365.25, type='', period=NULL, missRect=FALSE, time_header=NULL, info_header=TRUE, time_ratio=2, var_ratio=3) {
+panels_layout = function (df_data, df_meta, layout_matrix, figdir='', filedir_opt='', filename_opt='', variable='', df_trend=NULL, p_threshold=0.1, unit2day=365.25, type='', trend_period=NULL, mean_period=NULL, missRect=FALSE, time_header=NULL, info_header=TRUE, time_ratio=2, var_ratio=3) {
- if (all(class(df_data) != 'list')) {
- df_data = list(df_data)
+ outfile = "Panels"
+ if (filename_opt != '') {
+ outfile = paste(outfile, '_', filename_opt, sep='')
+ }
+ outfile = paste(outfile, '.pdf', sep='')
+
+ # If there is not a dedicated figure directory it creats one
+ outdir = file.path(figdir, filedir_opt, sep='')
+ if (!(file.exists(outdir))) {
+ dir.create(outdir)
+ }
+
+ outdirTmp = file.path(outdir, 'tmp')
+ if (!(file.exists(outdirTmp))) {
+ dir.create(outdirTmp)
+ } else {
+ unlink(outdirTmp, recursive=TRUE)
+ dir.create(outdirTmp)
}
nbp = length(df_data)
+ if (all(class(df_data) != 'list')) {
+ df_data = list(df_data)
+ }
+
if (all(class(df_trend) != 'list')) {
df_trend = list(df_trend)
if (length(df_trend) == 1) {
@@ -52,9 +72,37 @@ panels_layout = function (df_data, df_meta, layout_matrix, figdir='', filedir_op
missRect = replicate(nbp, missRect)
}}
+ # Get all different stations code
+ Code = levels(factor(df_meta$code))
+ nCode = length(Code)
+
+
+ # print(df_trend)
+ df_trendtmp = df_trend[[1]]
+
+ # print(df_trendtmp)
+
+ nPeriod_max = 0
+ for (code in Code) {
+
+ df_trend_code = df_trendtmp[df_trendtmp$code == code,]
+
+ 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))
+
+ if (nPeriod > nPeriod_max) {
+ nPeriod_max = nPeriod
+ }
+ }
+
list_df2plot = vector(mode='list', length=nbp)
- minTrend = c()
- maxTrend = c()
+ # minTrend = c()
+ # maxTrend = c()
for (i in 1:nbp) {
@@ -65,38 +113,92 @@ panels_layout = function (df_data, df_meta, layout_matrix, figdir='', filedir_op
type=type[[i]],
missRect=missRect[[i]])
- okTrend = df_trend[[i]]$trend[df_trend[[i]]$p <= p_threshold[[i]]]
+ # okTrend = df_trend[[i]]$trend[df_trend[[i]]$p <= p_threshold[[i]]]
- minTrend[i] = min(okTrend, na.rm=TRUE)
- maxTrend[i] = max(okTrend, na.rm=TRUE)
+ # minTrend[i] = min(okTrend, na.rm=TRUE)
+ # maxTrend[i] = max(okTrend, na.rm=TRUE)
list_df2plot[[i]] = df2plot
}
+ Start_code = vector(mode='list', length=nCode)
+ End_code = vector(mode='list', length=nCode)
+ Code_code = vector(mode='list', length=nCode)
+ Periods_code = vector(mode='list', length=nCode)
- outfile = "Panels"
- if (filename_opt != '') {
- outfile = paste(outfile, '_', filename_opt, sep='')
- }
- outfile = paste(outfile, '.pdf', sep='')
+ for (j in 1:nCode) {
+
+ code = Code[j]
- # If there is not a dedicated figure directory it creats one
- outdir = file.path(figdir, filedir_opt, sep='')
- if (!(file.exists(outdir))) {
- dir.create(outdir)
+ df_trend_code = df_trendtmp[df_trendtmp$code == code,]
+
+ 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 = c()
+
+ for (i in 1:nPeriod_max) {
+ Periods = append(Periods,
+ paste(substr(Start[i], 1, 4),
+ substr(End[i], 1, 4),
+ sep=' / '))
+ }
+ Start_code[[j]] = Start
+ End_code[[j]] = End
+ Code_code[[j]] = code
+ Periods_code[[j]] = Periods
}
- outdirTmp = file.path(outdir, 'tmp')
- if (!(file.exists(outdirTmp))) {
- dir.create(outdirTmp)
- } else {
- unlink(outdirTmp, recursive=TRUE)
- dir.create(outdirTmp)
+ TrendMean_code = array(rep(1, nPeriod_max*nbp*nCode),
+ dim=c(nPeriod_max, nbp, nCode))
+
+ for (j in 1:nPeriod_max) {
+
+ for (k in 1:nCode) {
+
+ code = Code[k]
+
+ for (i in 1:nbp) {
+
+ df_data = list_df2plot[[i]]$data
+ df_trend = list_df2plot[[i]]$trend
+ p_threshold = list_df2plot[[i]]$p_threshold
+
+ df_data_code = df_data[df_data$code == code,]
+ df_trend_code = df_trend[df_trend$code == code,]
+
+ Start = Start_code[Code_code == code][[1]][j]
+ End = End_code[Code_code == code][[1]][j]
+ Periods = Periods_code[Code_code == code][[1]][j]
+
+ df_data_code_per =
+ df_data_code[df_data_code$Date >= Start
+ & df_data_code$Date <= End,]
+
+ df_trend_code_per =
+ df_trend_code[df_trend_code$period_start == Start
+ & df_trend_code$period_end == End,]
+
+ Ntrend = nrow(df_trend_code_per)
+ if (Ntrend > 1) {
+ df_trend_code_per = df_trend_code_per[1,]
+ }
+
+ dataMean = mean(df_data_code_per$Qm3s, na.rm=TRUE)
+ trendMean = df_trend_code_per$trend / dataMean
+
+ TrendMean_code[j, i, k] = trendMean
+ }
+ }
}
- # Get all different stations code
- Code = levels(factor(df_meta$code))
- nCode = length(Code)
+ minTrendMean = apply(TrendMean_code, c(1, 2), min, na.rm=TRUE)
+ maxTrendMean = apply(TrendMean_code, c(1, 2), max, na.rm=TRUE)
for (code in Code) {
@@ -118,7 +220,7 @@ panels_layout = function (df_data, df_meta, layout_matrix, figdir='', filedir_op
time_header_code = time_header[time_header$code == code,]
Htime = time_panel(time_header_code, df_trend_code=NULL,
- period=period, missRect=TRUE,
+ trend_period=trend_period, missRect=TRUE,
unit2day=365.25, type='Q', first=FALSE)
P[[2]] = Htime
@@ -138,13 +240,41 @@ panels_layout = function (df_data, df_meta, layout_matrix, figdir='', filedir_op
df_trend_code = df_trend[df_trend$code == code,]
color = c()
- for (j in 1:nrow(df_trend_code)) {
+ # for (j in 1:nrow(df_trend_code)) {
+ for (j in 1:nPeriod_max) {
if (df_trend_code$p[j] <= p_threshold){
- color_res = get_color(df_trend_code$trend[j],
- minTrend[i],
- maxTrend[i],
+ # color_res = get_color(df_trend_code$trend[j],
+ # minTrend[i],
+ # maxTrend[i],
+ # palette_name='perso',
+ # reverse=TRUE)
+
+ Start = Start_code[Code_code == code][[1]][j]
+ End = End_code[Code_code == code][[1]][j]
+ Periods = Periods_code[Code_code == code][[1]][j]
+
+ df_data_code_per =
+ df_data_code[df_data_code$Date >= Start
+ & df_data_code$Date <= End,]
+
+ df_trend_code_per =
+ df_trend_code[df_trend_code$period_start == Start
+ & df_trend_code$period_end == End,]
+
+ Ntrend = nrow(df_trend_code_per)
+ if (Ntrend > 1) {
+ df_trend_code_per = df_trend_code_per[1,]
+ }
+
+ dataMean = mean(df_data_code$Qm3s, na.rm=TRUE)
+ trendMean = df_trend_code_per$trend / dataMean
+
+ color_res = get_color(trendMean,
+ minTrendMean[j, i],
+ maxTrendMean[j, i],
palette_name='perso',
reverse=TRUE)
+
colortmp = color_res$color
} else {
colortmp = NA
@@ -155,9 +285,9 @@ panels_layout = function (df_data, df_meta, layout_matrix, figdir='', filedir_op
p = time_panel(df_data_code, df_trend_code, type=type,
p_threshold=p_threshold, missRect=missRect,
- unit2day=unit2day, last=(i > nbp-nbcol),
- color=color)
-
+ mean_period=mean_period, unit2day=unit2day,
+ last=(i > nbp-nbcol), color=color)
+
P[[i+nbh]] = p
}
@@ -208,7 +338,7 @@ panels_layout = function (df_data, df_meta, layout_matrix, figdir='', filedir_op
}
# By
- slice = 12 ##
+ slice = 12
nMat = as.integer(nCode/slice) + 1
sublist_df2plot = list_df2plot
@@ -231,8 +361,8 @@ panels_layout = function (df_data, df_meta, layout_matrix, figdir='', filedir_op
}
mat = matrice_panel(sublist_df2plot,
- subdf_meta,
- period)
+ subdf_meta,
+ trend_period=trend_period)
# Saving matrix plot
ggsave(plot=mat,
diff --git a/plotting/panel.R b/plotting/panel.R
index f19ff0c3dd7aeef9fcd3cf579283fa74ee171078..22d332577e00bad2f149ad55f3ac0b2fa386e90b 100644
--- a/plotting/panel.R
+++ b/plotting/panel.R
@@ -10,7 +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, first=FALSE, color=NULL) {
+time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missRect=FALSE, unit2day=365.25, trend_period=NULL, mean_period=NULL, last=FALSE, first=FALSE, color=NULL) {
if (type == 'sqrt(Q)') {
df_data_code$Qm3s = sqrt(df_data_code$Qm3s)
@@ -56,16 +56,13 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR
theme(panel.background=element_rect(fill='white'),
text=element_text(family='sans'),
- # panel.border=element_blank(),
panel.border = element_rect(color="grey85",
fill=NA,
size=0.7),
- # panel.grid.major.y=element_line(color='grey85', size=0.3),
panel.grid.major.y=element_line(color='grey85', size=0.15),
panel.grid.major.x=element_blank(),
- # axis.ticks.y=element_blank(),
axis.ticks.y=element_line(color='grey75', size=0.3),
axis.ticks.x=element_line(color='grey75', size=0.3),
@@ -79,7 +76,7 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR
hjust=-1E-3, color='grey20'),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
- # axis.title.y=element_text(size=8, color='grey20'),
+
axis.line.x=element_blank(),
axis.line.y=element_blank(),
)
@@ -103,26 +100,26 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR
}
}
- if ((type == 'sqrt(Q)' | type == 'Q') & !is.null(period)) {
+ if ((type == 'sqrt(Q)' | type == 'Q') & !is.null(trend_period)) {
- period = as.list(period)
+ trend_period = as.list(trend_period)
Imin = 10^99
- for (per in period) {
+ for (per in trend_period) {
I = interval(per[1], per[2])
if (I < Imin) {
Imin = I
- period_min = as.Date(per)
+ trend_period_min = as.Date(per)
}
}
p = p +
geom_rect(aes(xmin=min(df_data_code$Date),
ymin=0,
- xmax=period_min[1],
+ xmax=trend_period_min[1],
ymax= maxQ*1.1),
linetype=0, fill='grey85', alpha=0.3) +
- geom_rect(aes(xmin=period_min[2],
+ geom_rect(aes(xmin=trend_period_min[2],
ymin=0,
xmax=max(df_data_code$Date),
ymax= maxQ*1.1),
@@ -137,7 +134,7 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR
} else {
p = p +
geom_point(aes(x=df_data_code$Date, y=df_data_code$Qm3s),
- shape=1, color='grey20', size=1)
+ shape=21, color='grey50', fill='grey95', size=1)
}
if (missRect) {
@@ -153,94 +150,114 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR
ymax=maxQ*1.1),
linetype=0, fill='Wheat', alpha=0.4)
}
+
+ if (!is.null(mean_period)) {
+ mean_period = as.list(mean_period)
+ nPeriod_mean = length(mean_period)
+
+ for (j in 1:nPeriod_mean) {
+ Start_mean = mean_period[[j]][1]
+ End_mean = mean_period[[j]][2]
+
+ df_data_code_per =
+ df_data_code[df_data_code$Date >= Start_mean
+ & df_data_code$Date <= End_mean,]
+
+ absMin = min(df_data_code_per$Date)
+ if (absMin > min(df_data_code$Date)) {
+ absMin = absMin - months(6)
+ }
+
+ absMax = max(df_data_code_per$Date)
+ if (absMax < max(df_data_code$Date)) {
+ absMax = absMax + months(6)
+ }
+
+ abs = c(absMin, absMax)
+ ord = rep(mean(df_data_code_per$Qm3s, na.rm=TRUE), times=2)
+ plot_mean = tibble(abs=abs, ord=ord)
+
+ p = p +
+ geom_line(data=plot_mean, aes(x=abs, y=ord),
+ color='grey40', linetype='solid',
+ size=0.5)
+ }
+ }
if (!is.null(df_trend_code)) {
- # print(df_trend_code)
-
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]))
- # }
+ nPeriod_trend = max(length(UStart), length(UEnd))
ltype = c('solid', 'dashed', 'dotted', 'twodash')
lty = c('solid', '22', 'dotted', 'twodash')
- ii = 0
- for (i in 1:nPeriod) {
+ # ii = 0
+ for (i in 1:nPeriod_trend) {
df_trend_code_per =
df_trend_code[df_trend_code$period_start == Start[i]
& df_trend_code$period_end == End[i],]
- if (df_trend_code_per$p <= p_threshold) {
-
- ii = ii + 1
-
- 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_num = as.numeric(abs) / unit2day
+ 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_num = as.numeric(abs) / unit2day
+ ord = abs_num * df_trend_code_per$trend +
+ df_trend_code_per$intercept
+ plot_trend = tibble(abs=abs, ord=ord)
+ if (df_trend_code_per$p <= p_threshold) {
+ p = p +
+ geom_line(data=plot_trend, aes(x=abs, y=ord),
+ color=color[i],
+ linetype=ltype[i], size=0.7)
+ colortxt = color[i]
+ } else {
+ p = p +
+ geom_line(data=plot_trend, aes(x=abs, y=ord),
+ color='grey85', linetype=ltype[i],
+ size=0.5)
+ colortxt = 'grey85'
+ }
- ord = abs_num * df_trend_code_per$trend +
- df_trend_code_per$intercept
+ codeDate = df_data_code$Date
+ codeQ = df_data_code$Qm3s
+
+ x = gpct(2, codeDate, shift=TRUE)
+ xend = x + gpct(3, codeDate)
+
+ dy = gpct(7, codeQ, ref=0)
+ y = gpct(100, codeQ, ref=0) - (i-1)*dy
- plot = tibble(abs=abs, ord=ord)
+ xt = xend + gpct(1, codeDate)
- if (!is.na(color[i])) {
- 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')
- }
+ trend = df_trend_code_per$trend
+ power = power = get_power(trend)
+ powerC = as.character(power)
+ brk = 10^power
+ trendC = as.character(round(trend / brk, 2))
- codeDate = df_data_code$Date
- codeQ = df_data_code$Qm3s
+ label = bquote(bold(.(trendC)~'x'~'10'^{.(powerC)})~'['*m^{3}*'.'*s^{-1}*'.'*an^{-1}*']')
+
+ p = p +
+ annotate("segment",
+ x=x, xend=xend,
+ y=y, yend=y,
+ color=colortxt,
+ lty=lty[i], lwd=1) +
- x = gpct(2, codeDate, shift=TRUE)
- xend = x + gpct(3, codeDate)
-
- dy = gpct(7, codeQ, ref=0)
- y = gpct(100, codeQ, ref=0) - (ii-1)*dy
-
- xt = xend + gpct(1, codeDate)
-
- trend = df_trend_code_per$trend
- power = power = get_power(trend)
- powerC = as.character(power)
- brk = 10^power
- trendC = as.character(round(trend / brk, 2))
-
- label = bquote(bold(.(trendC)~'x'~'10'^{.(powerC)})~'['*m^{3}*'.'*s^{-1}*'.'*an^{-1}*']')
-
- p = p +
- annotate("segment",
- x=x, xend=xend,
- y=y, yend=y,
- color=color[i],
- lty=lty[i], lwd=1) +
-
- annotate("text",
- label=label, size=3,
- x=xt, y=y,
- hjust=0, vjust=0.4,
- color=color[i])
- }
+ annotate("text",
+ label=label, size=3,
+ x=xt, y=y,
+ hjust=0, vjust=0.4,
+ color=colortxt)
}
}
@@ -362,7 +379,7 @@ text_panel = function(code, df_meta) {
-matrice_panel = function (list_df2plot, df_meta, period) {
+matrice_panel = function (list_df2plot, df_meta, trend_period) {
nbp = length(list_df2plot)
@@ -620,10 +637,10 @@ matrice_panel = function (list_df2plot, df_meta, period) {
x=x, xend=xend,
y=y, yend=yend,
color="grey40", size=0.35)
-
+
yt = y + 0.15
- Start = period[[j]][1]
- End = period[[j]][2]
+ Start = trend_period[[j]][1]
+ End = trend_period[[j]][2]
periodName = bquote(bold('Période')~bold(.(as.character(j))))
# bquote(bold(.(Start))~'/'~bold(.(End)))
@@ -757,7 +774,87 @@ matrice_panel = function (list_df2plot, df_meta, period) {
-histogram = function (data_bin, bins, binwidth=NULL, figdir='', filedir_opt='') {
+histogram = function (data_bin, df_meta, breaks=NULL, bins=NULL, binwidth=NULL, figdir='', filedir_opt='') {
+
+ # Get all different stations code
+ Code = levels(factor(df_meta$code))
+ nCode = length(Code)
+
+ # If there is not a dedicated figure directory it creats one
+ outdir = file.path(figdir, filedir_opt, sep='')
+ if (!(file.exists(outdir))) {
+ dir.create(outdir)
+ }
+
+ datebreak = 10
+ dateminbreak = 1
+
+ res_hist = hist(data_bin, breaks='years', plot=FALSE)
+ counts = res_hist$counts
+ counts_pct = counts/nCode * 100
+ breaks = as.Date(res_hist$breaks)
+ mids = as.Date(res_hist$mids)
+
+ p = ggplot() +
+
+ theme(panel.background=element_rect(fill='white'),
+ text=element_text(family='sans'),
+
+ panel.border = element_rect(color="grey85",
+ fill=NA,
+ size=0.7),
+
+ panel.grid.major.y=element_line(color='grey85', size=0.15),
+ panel.grid.major.x=element_blank(),
+
+ axis.ticks.y=element_line(color='grey75', size=0.3),
+ axis.ticks.x=element_line(color='grey75', size=0.3),
+
+ axis.text.x=element_text(color='grey40'),
+ axis.text.y=element_text(color='grey40'),
+
+ ggh4x.axis.ticks.length.minor=rel(0.5),
+ axis.ticks.length=unit(1.5, 'mm'),
+
+ plot.title=element_text(size=9, vjust=-2,
+ hjust=-1E-3, color='grey20'),
+ axis.title.x=element_blank(),
+ axis.title.y=element_blank(),
+ axis.line.x=element_blank(),
+ axis.line.y=element_blank(),
+ ) +
+
+ geom_bar(aes(x=mids, y=counts_pct),
+ stat='identity',
+ # color="#00A3A8",
+ fill="#00A3A8") +
+
+ scale_x_date(date_breaks=paste(as.character(datebreak),
+ 'year', sep=' '),
+ date_minor_breaks=paste(as.character(dateminbreak),
+ 'year', sep=' '),
+ guide='axis_minor',
+ date_labels="%Y",
+ limits=c(min(data_bin)-years(0),
+ max(data_bin)+years(0)),
+ expand=c(0, 0)) +
+
+ scale_y_continuous(limits=c(0,
+ max(counts_pct)*1.1),
+ expand=c(0, 0))
+
+ ggsave(plot=p,
+ path=outdir,
+ filename=paste('hist_break_date', '.pdf', sep=''),
+ width=10, height=10, units='cm', dpi=100)
+}
+
+
+cumulative = function (data_bin, df_meta, dyear=10, breaks=NULL, bins=NULL, binwidth=NULL, figdir='', filedir_opt='') {
+
+ # Get all different stations code
+ Code = levels(factor(df_meta$code))
+ nCode = length(Code)
# If there is not a dedicated figure directory it creats one
outdir = file.path(figdir, filedir_opt, sep='')
@@ -768,27 +865,45 @@ histogram = function (data_bin, bins, binwidth=NULL, figdir='', filedir_opt='')
datebreak = 10
dateminbreak = 1
- res = stat_bin(x=data_bin,
- bins=bins,
- binwidth=binwidth)
+ res_hist = hist(data_bin, breaks='years', plot=FALSE)
+ counts = res_hist$counts
+ cumul = cumsum(counts)
+ cumul_pct = cumul/nCode * 100
+ breaks = as.Date(res_hist$breaks)
+ mids = as.Date(res_hist$mids)
+
+ mids = c(mids[1] - years(dyear), mids[1] - years(1),
+ mids,
+ mids[length(mids)] + years(dyear))
+ cumul_pct = c(0, 0,
+ cumul_pct,
+ cumul_pct[length(cumul_pct)])
- print(res)
+ mids = mids + months(6)
+
+ breaks = breaks + 1
+ breaks = breaks[-length(breaks)]
+
+ DB = c()
+ for (i in 1:length(breaks)) {
+ DB = c(DB, rep(breaks[i], times=counts[i]))
+ }
+ q50 = as.Date(quantile(DB, probs=0.5)) + years(1)
+
+ print(paste('mediane :', q50))
p = ggplot() +
theme(panel.background=element_rect(fill='white'),
text=element_text(family='sans'),
- # panel.border=element_blank(),
panel.border = element_rect(color="grey85",
fill=NA,
size=0.7),
- # panel.grid.major.y=element_line(color='grey85', size=0.3),
panel.grid.major.y=element_line(color='grey85', size=0.15),
panel.grid.major.x=element_blank(),
- # axis.ticks.y=element_blank(),
axis.ticks.y=element_line(color='grey75', size=0.3),
axis.ticks.x=element_line(color='grey75', size=0.3),
@@ -802,37 +917,35 @@ histogram = function (data_bin, bins, binwidth=NULL, figdir='', filedir_opt='')
hjust=-1E-3, color='grey20'),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
- # axis.title.y=element_text(size=8, color='grey20'),
axis.line.x=element_blank(),
axis.line.y=element_blank(),
) +
- geom_histogram(aes(x=data_bin),
- bins=bins,
- binwidth=binwidth,
- color=NA,
- fill="#00A3A8") +
+ geom_line(aes(x=mids, y=cumul_pct),
+ color="#00A3A8") +
+ geom_line(aes(x=c(q50, q50), y=c(0, 100)),
+ color="wheat",
+ lty='dashed') +
+
scale_x_date(date_breaks=paste(as.character(datebreak),
'year', sep=' '),
date_minor_breaks=paste(as.character(dateminbreak),
'year', sep=' '),
guide='axis_minor',
date_labels="%Y",
- limits=c(min(data_bin),
- max(data_bin)),
- expand=c(0, 0)) #+
+ limits=c(min(mids)-years(0),
+ max(mids)+years(0)),
+ expand=c(0, 0)) +
- # scale_y_continuous(limits=c(0, ),
- # expand=c(0, 0))
+ scale_y_continuous(limits=c(-1, 101),
+ expand=c(0, 0))
ggsave(plot=p,
path=outdir,
- filename=paste('break_date', '.pdf', sep=''),
+ filename=paste('cumul_break_date', '.pdf', sep=''),
width=10, height=10, units='cm', dpi=100)
}
-
-
get_color = function (value, min, max, ncolor=256, palette_name='perso', reverse=FALSE) {
@@ -845,13 +958,25 @@ get_color = function (value, min, max, ncolor=256, palette_name='perso', reverse
# '#fdb147',
# '#fd4659'
- '#0f3b57',
- '#1d7881',
- '#80c4a9',
- '#e2dac6', #mid
- '#fadfad',
- '#d08363',
- '#7e392f'
+ # '#0f3b57',
+ # '#1d7881',
+ # '#80c4a9',
+ # '#e2dac6', #mid
+ # '#fadfad',
+ # '#d08363',
+ # '#7e392f'
+
+ '#193830',
+ '#2A6863',
+ '#449C93',
+ '#7ACEB9',
+ '#BCE6DB',
+ '#EFE0B0',
+ '#D4B86A',
+ '#B3762A',
+ '#7F4A23',
+ '#452C1A'
+
))(ncolor)
} else {
@@ -897,7 +1022,6 @@ void = ggplot() + geom_blank(aes(1,1)) +
)
-
palette_tester = function () {
n = 300
@@ -923,14 +1047,25 @@ palette_tester = function () {
# '#01665e',
# '#003c30'
- '#0f3b57',
- '#1d7881',
- '#80c4a9',
- '#e2dac6', #mid
- '#fadfad',
- '#d08363',
- '#7e392f'
+ # '#0f3b57',
+ # '#1d7881',
+ # '#80c4a9',
+ # '#e2dac6', #mid
+ # '#fadfad',
+ # '#d08363',
+ # '#7e392f'
+ '#193830',
+ '#2A6863',
+ '#449C93',
+ '#7ACEB9',
+ '#BCE6DB',
+ '#EFE0B0',
+ '#D4B86A',
+ '#B3762A',
+ '#7F4A23',
+ '#452C1A'
+
))(n)
p = ggplot() +
diff --git a/processing/analyse.R b/processing/analyse.R
index 85a9f67ae3c7062de5519e8766a36494ca4927b7..d86cb80086eaadda821c2b96fe6dd36776f4f642 100644
--- a/processing/analyse.R
+++ b/processing/analyse.R
@@ -142,31 +142,37 @@ get_period = function (per, df_Xtrend, df_XEx, df_Xlist) {
-get_break = function (df_data, df_meta) {
+get_break = function (df_data, df_meta, p_thresold=0.05) {
# Get all different stations code
Code = levels(factor(df_meta$code))
nCode = length(Code)
date_break = list()
+ Code_break = c()
for (code in Code) {
df_data_code = df_data[df_data$code == code,]
df_data_codeNoNA = df_data_code[!is.na(df_data_code$Qm3s),]
res_break = pettitt.test(df_data_codeNoNA$Qm3s)
+
+ p_value = res_break$p
nbreak = res_break$nobs
ibreak = res_break$estimate
if (length(ibreak) > 1) {
ibreak = ibreak[1]
}
- step1 = mean(df_data_codeNoNA$Qm3s[1:ibreak])
- step2 = mean(df_data_codeNoNA$Qm3s[(ibreak+1):nbreak])
- date_break = append(date_break, df_data_codeNoNA$Date[ibreak])
+ # step1 = mean(df_data_codeNoNA$Qm3s[1:ibreak])
+ # step2 = mean(df_data_codeNoNA$Qm3s[(ibreak+1):nbreak])
+ if (p_value <= p_thresold) {
+ date_break = append(date_break,
+ df_data_codeNoNA$Date[ibreak])
+ Code_break = append(Code_break, code)
+ }
}
-
- df_break = tibble(code=Code, Date=as.Date(date_break))
+ df_break = tibble(code=Code_break, Date=as.Date(date_break))
return (df_break)
}
diff --git a/script.R b/script.R
index e20f9331b77bcb35fcdeba912ad69fcfba40bdc2..a8c49488acb93fa1d75c5b75da5793b7f5a8c428 100644
--- a/script.R
+++ b/script.R
@@ -23,7 +23,7 @@ filedir =
### MANUAL SELECTION ###
# Name of the file that will be analysed from the AG directory
filename =
- ""
+ # ""
# c(
# "S2235610_HYDRO_QJM.txt",
@@ -33,9 +33,9 @@ filename =
# "A2250310_HYDRO_QJM.txt"
# )
- # c("O3035210_HYDRO_QJM.txt",
- # "O3011010_HYDRO_QJM.txt",
- # "O1442910_HYDRO_QJM.txt")
+ c("O3035210_HYDRO_QJM.txt",
+ "O3011010_HYDRO_QJM.txt",
+ "O1442910_HYDRO_QJM.txt")
@@ -45,8 +45,8 @@ AGlistdir =
""
AGlistname =
- # ""
- "Liste-station_RRSE.docx"
+ ""
+ # "Liste-station_RRSE.docx"
### NIVALE SELECTION ###
@@ -61,10 +61,14 @@ NVlistname =
### TREND ANALYSIS ###
# Time period to analyse
-period_all = c("1800-01-01", "2019-12-31")
-period1 = c("1968-01-01", "1988-12-31")
-period2 = c("1988-01-01", "2019-12-31")
-period = list(period_all, period1, period2)
+periodAll = c("1800-01-01", "2019-12-31")
+periodSub = c("1968-01-01", "2019-12-31")
+trend_period = list(periodAll, periodSub)
+
+# Time period to mean
+period1 = c("1968-01-01", "1994-12-31")
+period2 = c("1995-01-01", "2019-12-31")
+mean_period = list(period1, period2)
# p value
p_thresold = 0.1 #c(0.01, 0.05, 0.1)
@@ -184,22 +188,32 @@ df_meta = df_join$meta
# QA TREND #
-res_QAtrend = get_QAtrend(df_data, period=period, p_thresold=p_thresold)
+res_QAtrend = get_QAtrend(df_data, period=trend_period,
+ p_thresold=p_thresold)
# QMNA TREND #
-# res_QMNAtrend = get_QMNAtrend(df_data, period=period,
-# p_thresold=p_thresold)
+res_QMNAtrend = get_QMNAtrend(df_data, period=trend_period,
+ p_thresold=p_thresold)
+
+# VCN10 TREND #
+res_VCN10trend = get_VCN10trend(df_data, df_meta,
+ period=trend_period,
+ p_thresold=p_thresold)
-# # VCN10 TREND #
-# res_VCN10trend = get_VCN10trend(df_data, df_meta,
-# period=period, p_thresold=p_thresold)
+# df_break = get_break(res_QAtrend$data, df_meta)
+# df_break = get_break(res_QMNAtrend$data, df_meta)
+# df_break = get_break(res_VCN10trend$data, df_meta)
-df_break = get_break(res_QAtrend$data, df_meta)
-histogram(df_break$Date,
- bins=nrow(df_break),
- # binwidth=duration(1, units='year'),
- figdir=figdir)
+# histogram(df_break$Date, df_meta,
+ # breaks=seq(min(df_break$Date),
+ # max(df_break$Date), "years"),
+ # figdir=figdir)
+
+# cumulative(df_break$Date, df_meta, dyear=8,
+ # breaks=seq(min(df_break$Date),
+ # max(df_break$Date), "years"),
+ # figdir=figdir)
# TIME PANEL #
# Plot time panel of debit by stations
@@ -214,22 +228,22 @@ histogram(df_break$Date,
# figdir=figdir,
# filename_opt='time')
-# panels_layout(list(res_QAtrend$data, res_QMNAtrend$data,
-# res_VCN10trend$data),
-# layout_matrix=c(1, 2, 3),
-# df_meta=df_meta,
-# df_trend=list(res_QAtrend$trend, res_QMNAtrend$trend,
-# res_VCN10trend$trend),
-# type=list(bquote(Q[A]), bquote(Q[MNA]), bquote(V[CN10])),
-# missRect=list(TRUE, TRUE, TRUE),
-# period=period,
-# info_header=TRUE,
-# time_header=df_data,
-# time_ratio=2,
-# var_ratio=3,
-# figdir=figdir,
-# filename_opt='')
-
+panels_layout(list(res_QAtrend$data, res_QMNAtrend$data,
+ res_VCN10trend$data),
+ layout_matrix=c(1, 2, 3),
+ df_meta=df_meta,
+ df_trend=list(res_QAtrend$trend, res_QMNAtrend$trend,
+ res_VCN10trend$trend),
+ type=list(bquote(Q[A]), bquote(Q[MNA]), bquote(V[CN10])),
+ missRect=list(TRUE, TRUE, TRUE),
+ trend_period=trend_period,
+ mean_period=mean_period,
+ info_header=TRUE,
+ time_header=df_data,
+ time_ratio=2,
+ var_ratio=3,
+ figdir=figdir,
+ filename_opt='')
# panels_layout(list(res_QAtrend$data, res_VCN10trend$data),
# layout_matrix=c(1, 2),
@@ -238,7 +252,7 @@ histogram(df_break$Date,
# res_VCN10trend$trend),
# type=list(bquote(Q[A]), bquote(V[CN10])),
# missRect=list(TRUE, TRUE),
-# period=period,
+# period=trend_period,
# info_header=TRUE,
# time_header=df_data,
# time_ratio=2,
@@ -254,7 +268,7 @@ histogram(df_break$Date,
# df_trend=list(res_QMNAtrend$trend),
# type=list(bquote(Q[MNA])),
# missRect=list(TRUE),
-# period=period,
+# period=trend_period,
# info_header=TRUE,
# time_header=df_data,
# time_ratio=2,