diff --git a/plotting/layout.R b/plotting/layout.R index fcb7b4a2a7a54913089f954644ad7d47f2b25c97..f5072a77df3b5230ccd0ea1182ac989dbb40face 100644 --- a/plotting/layout.R +++ b/plotting/layout.R @@ -144,7 +144,7 @@ panels_layout = function (df_data, df_meta, layout_matrix, figdir='', filedir_op minTrend[i], maxTrend[i], palette_name='perso', - reverse=FALSE) + reverse=TRUE) colortmp = color_res$color } else { colortmp = NA diff --git a/plotting/panel.R b/plotting/panel.R index 88c470a37654c0bd642604b5e82eab294547cde6..f19ff0c3dd7aeef9fcd3cf579283fa74ee171078 100644 --- a/plotting/panel.R +++ b/plotting/panel.R @@ -101,33 +101,7 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR p = p + theme(plot.margin=margin(0, 5, 0, 5, unit="mm")) } - } - - - if (type == 'sqrt(Q)' | type == 'Q') { - p = p + - 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$Qm3s), - shape=1, color='grey20', size=1) - } - - if (missRect) { - 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] - - p = p + - geom_rect(aes(xmin=NAdate_Down, - ymin=0, - xmax=NAdate_Up, - ymax=maxQ*1.1), - linetype=0, fill='Wheat', alpha=0.4) - } + } if ((type == 'sqrt(Q)' | type == 'Q') & !is.null(period)) { @@ -155,6 +129,30 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR linetype=0, fill='grey85', alpha=0.3) } + if (type == 'sqrt(Q)' | type == 'Q') { + p = p + + 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$Qm3s), + shape=1, color='grey20', size=1) + } + + if (missRect) { + 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] + + p = p + + geom_rect(aes(xmin=NAdate_Down, + ymin=0, + xmax=NAdate_Up, + ymax=maxQ*1.1), + linetype=0, fill='Wheat', alpha=0.4) + } if (!is.null(df_trend_code)) { @@ -217,7 +215,7 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR x = gpct(2, codeDate, shift=TRUE) xend = x + gpct(3, codeDate) - dy = gpct(6, codeQ, ref=0) + dy = gpct(7, codeQ, ref=0) y = gpct(100, codeQ, ref=0) - (ii-1)*dy xt = xend + gpct(1, codeDate) @@ -519,9 +517,10 @@ matrice_panel = function (list_df2plot, df_meta, period) { if (df_trend_code_per$p <= p_threshold){ color_res = get_color(trendMean, minTrendMean[j, i], - minTrendMean[j, i], + maxTrendMean[j, i], palette_name='perso', - reverse=FALSE) + reverse=TRUE) + fill = color_res$color color = 'white' Pthresold = p_thresold @@ -600,6 +599,9 @@ matrice_panel = function (list_df2plot, df_meta, period) { Fill_mat_per = Fill_mat[NPeriod_mat == j] Color_mat_per = Color_mat[NPeriod_mat == j] + # print(j) + # print(Fill_mat_per) + Xtmp = as.integer(factor(as.character(Type_mat_per))) Xc = j + (j - 1)*nbp*2 @@ -669,7 +671,7 @@ matrice_panel = function (list_df2plot, df_meta, period) { } mat = mat + - annotate('text', x=Xc, y=max(Y) + 0.8, + annotate('text', x=Xc, y=max(Y) + 0.85, label=bquote(bold('Début')), hjust=0.5, vjust=0.5, size=3, color='grey20') + @@ -702,27 +704,21 @@ matrice_panel = function (list_df2plot, df_meta, period) { if (nchar(name) > ncharMax) { name = paste(substr(name, 1, ncharMax), '...', sep='') } - label = bquote(.(name)~'-'~bold(.(code))) + # label = bquote(.(name)~'-'~bold(.(code))) # period_code = Periods_code[[i]] # nPeriod_code = length(period_code) - # dy = 0.17 - # shift = as.integer((nPeriod_code + 1) / 2) * dy ### not ok mat = mat + - annotate('text', x=0.3, y=i, - label=label, + annotate('text', x=0.3, y=i + 0.14, + label=bquote(bold(.(code))), + hjust=1, vjust=0.5, + size=3.5, color="#00A3A8") + + + annotate('text', x=0.3, y=i - 0.14, + label=name, hjust=1, vjust=0.5, size=3.5, color="#00A3A8") - # for (j in 1:nPeriod_code) { - # periodj = bquote('Période'~.(j)~':'~.(period_code[j])) - - # mat = mat + - # annotate('text', x=0.3, y=i + shift - dy*j, - # label=periodj, - # hjust=1, vjust=0.5, - # size=3, color="#00A3A8") - # } for (j in 1:nPeriod_max) { Xc = j + (j - 1)*nbp*2 @@ -731,12 +727,12 @@ matrice_panel = function (list_df2plot, df_meta, period) { periodEnd = substr(label, 8, 11) mat = mat + - annotate('text', x=Xc, y=i + 0.1, + annotate('text', x=Xc, y=i + 0.13, label=bquote(bold(.(periodStart))), hjust=0.5, vjust=0.5, size=3, color='grey40') + - annotate('text', x=Xc, y=i - 0.1, + annotate('text', x=Xc, y=i - 0.13, label=bquote(bold(.(periodEnd))), hjust=0.5, vjust=0.5, size=3, color='grey40') @@ -753,7 +749,7 @@ matrice_panel = function (list_df2plot, df_meta, period) { expand=c(0, 0)) + scale_y_continuous(limits=c(1 - rel(0.5), - height + rel(1.25)), + height + rel(1.5)), expand=c(0, 0)) return (mat) @@ -761,17 +757,101 @@ matrice_panel = function (list_df2plot, df_meta, period) { +histogram = function (data_bin, bins, binwidth=NULL, figdir='', filedir_opt='') { + + # 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 = stat_bin(x=data_bin, + bins=bins, + binwidth=binwidth) + + print(res) + + 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), + + 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.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") + + + 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)) #+ + + # scale_y_continuous(limits=c(0, ), + # expand=c(0, 0)) + + ggsave(plot=p, + path=outdir, + filename=paste('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) { if (palette_name == 'perso') { palette = colorRampPalette(c( - '#1a4157', - '#00af9d', - '#fbdd7e', - '#fdb147', - '#fd4659' + # '#1a4157', + # '#00af9d', + # '#fbdd7e', + # '#fdb147', + # '#fd4659' + + '#0f3b57', + '#1d7881', + '#80c4a9', + '#e2dac6', #mid + '#fadfad', + '#d08363', + '#7e392f' ))(ncolor) } else { @@ -825,24 +905,59 @@ palette_tester = function () { Y = rep(0, times=n) palette = colorRampPalette(c( - '#1a4157', - '#00af9d', - '#fbdd7e', - '#fdb147', - '#fd4659' + # '#1a4157', + # '#00af9d', + # '#fbdd7e', + # '#fdb147', + # '#fd4659' + + # '#543005', + # '#8c510a', + # '#bf812d', + # '#dfc27d', + # '#f6e8c3', + # '#f5f5f5', + # '#c7eae5', + # '#80cdc1', + # '#35978f', + # '#01665e', + # '#003c30' + + '#0f3b57', + '#1d7881', + '#80c4a9', + '#e2dac6', #mid + '#fadfad', + '#d08363', + '#7e392f' + ))(n) p = ggplot() + - geom_line(aes(x=X, y=Y), color=palette[X], size=10) + + + theme( + plot.background = element_blank(), + panel.grid.major = element_blank(), + panel.grid.minor = element_blank(), + panel.border = element_blank(), + panel.background = element_blank(), + axis.title.x = element_blank(), + axis.title.y = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + axis.line = element_blank() + ) + + + geom_line(aes(x=X, y=Y), color=palette[X], size=60) + scale_y_continuous(expand=c(0, 0)) ggsave(plot=p, - path='/figures', filename=paste('palette_test', '.pdf', sep=''), width=10, height=10, units='cm', dpi=100) } -# palette_teste() +# palette_tester() get_power = function (value) { diff --git a/processing/analyse.R b/processing/analyse.R index 8e2453dd7c89145993c1350e0d0bc9172c293c97..85a9f67ae3c7062de5519e8766a36494ca4927b7 100644 --- a/processing/analyse.R +++ b/processing/analyse.R @@ -3,6 +3,7 @@ library(dplyr) library(zoo) library(StatsAnalysisTrend) library(lubridate) +library(trend) # Sourcing R file @@ -140,49 +141,35 @@ get_period = function (per, df_Xtrend, df_XEx, df_Xlist) { } -# get_QAtrend = function (df_data, period, p_thresold) { -# # AVERAGE ANNUAL FLOW : QA # - -# period = as.list(period) -# for (p_thr in p_thresold) { +get_break = function (df_data, df_meta) { + + # Get all different stations code + Code = levels(factor(df_meta$code)) + nCode = length(Code) -# Imax = 0 -# df_QAtrendB = tibble() + date_break = list() + for (code in Code) { -# 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) - -# } + 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) + nbreak = res_break$nobs + ibreak = res_break$estimate -# } - -# res_QAtrend = clean(df_QAtrendB, df_QAExB, df_QAlistB) + 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]) + } + + df_break = tibble(code=Code, Date=as.Date(date_break)) -# return (res_QAtrend) -# } + return (df_break) +} diff --git a/processing/extract.R b/processing/extract.R index b7835335745c4b9efde67ff754a98e8956e40f88..e237006a8eed10d600caf40366db56a21c1f10c6 100644 --- a/processing/extract.R +++ b/processing/extract.R @@ -84,14 +84,22 @@ iRegHydro = c('D'='Affluents du Rhin', create_selection = function (computer_data_path, filedir, outname) { - + + # Out file for store results outfile = file.path(computer_data_path, outname) - codelist = c() - dir_path = file.path(computer_data_path, filedir) + + # Path to find the directory of desired codes + dir_path = file.path(computer_data_path, filedir) + # Create a filelist of all the filename in the above directory filelist_tmp = list.files(dir_path) - + # Create a filelist to store all station codes + codelist = c() + + # For all the filename in the file list for (f in filelist_tmp) { + # If the filename is a 'txt' file if (file_ext(f) == 'txt') { + # Then the station code is stored codelist = c(codelist, gsub('.txt', '', f)) } } diff --git a/script.R b/script.R index 1263530e7bba5ff8e063b798ba1f81ba599e9fdd..e20f9331b77bcb35fcdeba912ad69fcfba40bdc2 100644 --- a/script.R +++ b/script.R @@ -7,14 +7,14 @@ computer_data_path = "/home/louis/Documents/bouleau/INRAE/CDD_stationnarite/data" # "C:\\Users\\louis.heraut\\Documents\\CDD_stationnarite\\data" -# Work path +# Work path (it needs to end with '/ASH' directory) computer_work_path = "/home/louis/Documents/bouleau/INRAE/CDD_stationnarite/ASH" # "C:\\Users\\louis.heraut\\Documents\\CDD_stationnarite\\ASH" ### BANQUE HYDRO ### -# Path to the directory where BH data is stored +# Path to the directory where BH data is stored from the work path filedir = # "" "BanqueHydro_Export2021" @@ -39,7 +39,7 @@ filename = -### AGENCE ADOUR GARONNE SELECTION ### +### AGENCE EAU ADOUR GARONNE SELECTION ### # Path to the list file of AG data that will be analysed AGlistdir = "" @@ -62,8 +62,9 @@ NVlistname = ### TREND ANALYSIS ### # Time period to analyse period_all = c("1800-01-01", "2019-12-31") -period2 = c("1968-01-01", "2019-12-31") -period = list(period_all, period2) +period1 = c("1968-01-01", "1988-12-31") +period2 = c("1988-01-01", "2019-12-31") +period = list(period_all, period1, period2) # p value p_thresold = 0.1 #c(0.01, 0.05, 0.1) @@ -83,8 +84,6 @@ source('processing/analyse.R', encoding='latin1') source('plotting/panel.R', encoding='latin1') source('plotting/layout.R', encoding='latin1') -# Usefull library - # Result directory resdir = file.path(computer_work_path, 'results') @@ -101,13 +100,14 @@ if (!(file.exists(figdir))) { print(paste('figdir :', figdir)) +# Initialization of null data frame if there is no data selected df_data_AG = NULL df_data_NV = NULL df_meta_AG = NULL df_meta_NV = NULL -# AGENCE ADOUR GARONNE SELECTION # +# AGENCE EAU ADOUR GARONNE SELECTION # if (AGlistname != ""){ # Get only the selected station from a list station file @@ -123,15 +123,14 @@ if (AGlistname != ""){ 'choix'), c_num=c('BV_km2', 'longueur_serie')) - + + # Get filenames of the selection filename = df_selec_AG[df_selec_AG$ok,]$filename - ##### # filename = filename[(1+30):(16+30)] ##### - # Extract metadata about selected stations df_meta_AG = extract_meta(computer_data_path, filedir, filename) @@ -147,14 +146,13 @@ if (NVlistname != ""){ NVlistdir, NVlistname) + # Get filenames of the selection filename = df_selec_NV[df_selec_NV$ok,]$filename - ##### # filename = filename[(1+20):(16+20)] ##### - # Extract metadata about selected stations df_meta_NV = extract_meta(computer_data_path, filedir, filename) @@ -174,6 +172,7 @@ if (AGlistname == "" & NVlistname == "") { # JOIN # +# Make the join df_join = join(df_data_AG, df_data_NV, df_meta_AG, df_meta_NV) df_data = df_join$data df_meta = df_join$meta @@ -188,13 +187,19 @@ df_meta = df_join$meta res_QAtrend = get_QAtrend(df_data, period=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=period, +# p_thresold=p_thresold) + +# # VCN10 TREND # +# res_VCN10trend = get_VCN10trend(df_data, df_meta, +# period=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) +histogram(df_break$Date, + bins=nrow(df_break), + # binwidth=duration(1, units='year'), + figdir=figdir) # TIME PANEL # # Plot time panel of debit by stations @@ -209,21 +214,21 @@ res_VCN10trend = get_VCN10trend(df_data, df_meta, # 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), +# 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_VCN10trend$data), @@ -253,8 +258,8 @@ panels_layout(list(res_QAtrend$data, res_QMNAtrend$data, # info_header=TRUE, # time_header=df_data, # time_ratio=2, -# var_ratio=5, -# figdir=figdir, + # var_ratio=5, + # figdir=figdir, # filename_opt='') diff --git a/script_install.R b/script_install.R index 055e3dc44dc7ad57df5e0383283d24f4ae3bab94..f5b6b55f2b0b37749bb764d9b224f21fb82b5a23 100644 --- a/script_install.R +++ b/script_install.R @@ -13,6 +13,8 @@ install.packages("lubridate") install.packages('ggh4x') install.packages('extrafont') install.packages("RColorBrewer") +install.packages('trend') + library(devtools) install_github("https://github.com/benRenard/BFunk") #type '1'