From e2ce457b4bb5ecbe4fde5c845584e32466c5c9a0 Mon Sep 17 00:00:00 2001 From: "louis.heraut" <louis.heraut@inrae.fr> Date: Thu, 25 Nov 2021 18:59:02 +0100 Subject: [PATCH] Plot aes matrix --- plotting/panel.R | 433 ++++++++++++++++------------------------- processing/extractBH.R | 36 ++++ processing/extractNV.R | 38 +++- script.R | 50 ++--- 4 files changed, 266 insertions(+), 291 deletions(-) diff --git a/plotting/panel.R b/plotting/panel.R index d7d458e..49642ea 100644 --- a/plotting/panel.R +++ b/plotting/panel.R @@ -10,216 +10,7 @@ library(ggh4x) library(RColorBrewer) -# Time panel -panel = 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, header_ratio=2) { - - if (all(class(df_data) != 'list')) { - df_data = list(df_data) - } - - nbp = length(df_data) - - if (all(class(df_trend) != 'list')) { - df_trend = list(df_trend) - if (length(df_trend) == 1) { - df_trend = replicate(nbp, df_trend) - }} - - if (all(class(p_threshold) != 'list')) { - p_threshold = list(p_threshold) - if (length(p_threshold) == 1) { - p_threshold = replicate(nbp, p_threshold) - }} - - if (all(class(unit2day) != 'list')) { - unit2day = list(unit2day) - if (length(unit2day) == 1) { - unit2day = replicate(nbp, unit2day) - }} - - if (all(class(type) != 'list')) { - type = list(type) - if (length(type) == 1) { - type = replicate(nbp, type) - }} - - if (all(class(missRect) != 'list')) { - missRect = list(missRect) - if (length(missRect) == 1) { - missRect = replicate(nbp, missRect) - }} - - list_df2plot = vector(mode='list', length=nbp) - minTrend = c() - maxTrend = c() - nokTrend = c() - - for (i in 1:nbp) { - - df2plot = list(data=df_data[[i]], - trend=df_trend[[i]], - p_threshold=p_threshold[[i]], - unit2day=unit2day[[i]], - type=type[[i]], - missRect=missRect[[i]]) - - okTrend = df_trend[[i]]$p[df_trend[[i]]$p <= p_threshold[[i]]] - - print(okTrend) - - minTrend[i] = min(okTrend, na.rm=TRUE) - maxTrend[i] = max(okTrend, na.rm=TRUE) - nokTrend[i] = length(okTrend) - - list_df2plot[[i]] = df2plot - } - - - 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) - } - - - # Get all different stations code - Code = levels(factor(df_meta$code)) - nCode = length(Code) - - for (code in Code) { - - # Print code of the station for the current plotting - print(paste("Plotting for sation :", code)) - - nbh = as.numeric(info_header) + as.numeric(!is.null(time_header)) - nbg = nbp + nbh - - P = vector(mode='list', length=nbg) - - if (info_header) { - Htext = text_panel(code, df_meta) - P[[1]] = Htext - } - - if (!is.null(time_header)) { - time_header_code = time_header[time_header$code == code,] - - Htime = time_panel(time_header_code, df_trend_code=NULL, - period=period, missRect=TRUE, - unit2day=365.25, type='Q') - - P[[2]] = Htime - } - - - nbcol = ncol(as.matrix(layout_matrix)) - for (i in 1:nbp) { - df_data = list_df2plot[[i]]$data - df_trend = list_df2plot[[i]]$trend - p_threshold = list_df2plot[[i]]$p_threshold - unit2day = list_df2plot[[i]]$unit2day - missRect = list_df2plot[[i]]$missRect - type = list_df2plot[[i]]$type - - df_data_code = df_data[df_data$code == code,] - df_trend_code = df_trend[df_trend$code == code,] - - if (df_trend_code$p <= p_threshold){ - color_res = get_color(df_trend_code$p, - minTrend[i], - maxTrend[i], - ncolor=10, - palette_name="RdYlBu", - reverse=TRUE) - - color = color_res$color - palette = color_res$palette - - } else { - color = NULL - palette = NULL - } - - print(paste('min', minTrend[i])) - print(df_trend_code$p) - print(paste('max', maxTrend[i])) - - if (i == 1) {print(palette)} - - print(paste('color', color)) - print("") - - 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) - - P[[i+nbh]] = p - } - - layout_matrix = as.matrix(layout_matrix) - nel = nrow(layout_matrix)*ncol(layout_matrix) - - ## - idNA = which(is.na(layout_matrix), arr.ind=TRUE) - - layout_matrix[idNA] = seq(max(layout_matrix, na.rm=TRUE) + 1, - max(layout_matrix, na.rm=TRUE) + 1 + - nel) - ## - - layout_matrix_H = layout_matrix + nbh - - - LM = c() - LMcol = ncol(layout_matrix_H) - LMrow = nrow(layout_matrix_H) - for (i in 1:(LMrow+nbh)) { - - if (i <= nbh) { - LM = rbind(LM, rep(i, times=LMcol)) - } else { - LM = rbind(LM, - matrix(rep(layout_matrix_H[i-nbh,], - times=header_ratio), - ncol=LMcol, byrow=TRUE)) - }} - - plot = grid.arrange(grobs=P, layout_matrix=LM) - - # plot = grid.arrange(rbind(cbind(ggplotGrob(P[[2]]), ggplotGrob(P[[2]])), cbind(ggplotGrob(P[[3]]), ggplotGrob(P[[3]]))), heights=c(1/3, 2/3)) - - - # Saving - ggsave(plot=plot, - path=outdirTmp, - filename=paste(as.character(code), '.pdf', sep=''), - width=21, height=29.7, units='cm', dpi=100) - - } - - pdf_combine(input=file.path(outdirTmp, list.files(outdirTmp)), - output=file.path(outdir, outfile)) - unlink(outdirTmp, recursive=TRUE) -} - - - - - - -time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missRect=FALSE, unit2day=365.25, period=NULL, norm=TRUE, last=FALSE, color=NULL) { +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)') { @@ -238,33 +29,31 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR } dbrk = 10^power - ### /!\ PROBLÈME entre 11 et 10 ### - if (norm) { - df_data_code$Qm3s = df_data_code$Qm3s / dbrk - - if (!is.null(df_trend_code)) { - df_trend_code$p = df_trend_code$p / dbrk - df_trend_code$intercept = df_trend_code$intercept / dbrk - } + df_data_code$Qm3sN = df_data_code$Qm3s / dbrk - maxQ = max(df_data_code$Qm3s, na.rm=TRUE) + if (!is.null(df_trend_code)) { - if (maxQ >= 5) { - dbrk = 1.0 - accuracy = 0.1 - } else if (maxQ < 5 & maxQ >= 3) { - dbrk = 0.5 - accuracy = 0.1 - } else if (maxQ < 3 & maxQ >= 2) { - dbrk = 0.4 - accuracy = 0.1 - } else if (maxQ < 2 & maxQ >= 1) { - dbrk = 0.2 - accuracy = 0.1 - } else if (maxQ < 1) { - dbrk = 0.1 - accuracy = 0.1 - } + 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) + + if (maxQN >= 5) { + dbrk = 1.0 + accuracy = 0.1 + } else if (maxQN < 5 & maxQN >= 3) { + dbrk = 0.5 + accuracy = 0.1 + } else if (maxQN < 3 & maxQN >= 2) { + dbrk = 0.4 + accuracy = 0.1 + } else if (maxQN < 2 & maxQN >= 1) { + dbrk = 0.2 + accuracy = 0.1 + } else if (maxQN < 1) { + dbrk = 0.1 + accuracy = 0.1 } dDate = as.numeric(df_data_code$Date[length(df_data_code$Date)] - @@ -286,8 +75,6 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR p = ggplot() + # theme_bw() + - - ggtitle(bquote(.(type)~' ['*m^{3}*'.'*s^{-1}*'] x'~10^{.(as.character(power))})) + theme(panel.background=element_rect(fill='white'), text=element_text(family='sans'), @@ -325,19 +112,19 @@ 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$Qm3s), + geom_line(aes(x=df_data_code$Date, y=df_data_code$Qm3sN), color='grey20', size=0.3) } else { p = p + - # geom_line(aes(x=df_data_code$Date, y=df_data_code$Qm3s), + # geom_line(aes(x=df_data_code$Date, y=df_data_code$Qm3sN), # color='grey70') + - geom_point(aes(x=df_data_code$Date, y=df_data_code$Qm3s), + geom_point(aes(x=df_data_code$Date, y=df_data_code$Qm3sN), shape=1, color='grey20', size=1) } if (missRect) { - NAdate = df_data_code$Date[is.na(df_data_code$Qm3s)] + NAdate = df_data_code$Date[is.na(df_data_code$Qm3sN)] dNAdate = diff(NAdate) NAdate_Down = NAdate[append(Inf, dNAdate) != 1] NAdate_Up = NAdate[append(dNAdate, Inf) != 1] @@ -346,8 +133,8 @@ 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=maxQ*1.1), - linetype=0, fill='Wheat', alpha=0.3) + ymax=maxQN*1.1), + linetype=0, fill='Wheat', alpha=0.4) } if ((type == 'sqrt(Q)' | type == 'Q') & !is.null(period)) { @@ -356,13 +143,13 @@ 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= maxQ*1.1), + ymax= maxQN*1.1), linetype=0, fill='grey85', alpha=0.3) + geom_rect(aes(xmin=period[2], ymin=0, xmax=max(df_data_code$Date), - ymax= maxQ*1.1), + ymax= maxQN*1.1), linetype=0, fill='grey85', alpha=0.3) } @@ -374,19 +161,31 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR abs_num = as.numeric(abs) / unit2day - ord = abs_num * df_trend_code$trend + - df_trend_code$intercept + ord = abs_num * df_trend_code$trendN + + df_trend_code$interceptN if (!is.null(color)) { p = p + geom_line(aes(x=abs, y=ord), - color=color) + color=color, + size=0.7) } else { p = p + geom_line(aes(x=abs, y=ord), color='cornflowerblue') } - }} + + p = p + + ggtitle(bquote(.(type)~~'['*m^{3}*'.'*s^{-1}*'] x'~10^{.(as.character(power))}~~~'tendance :'~.(format(df_trend_code$trend, scientific=TRUE, digits=3))~m^{3}*'.'*s^{-1}*'.'*an^{-1})) + + } else { + p = p + + ggtitle(bquote(.(type)~' ['*m^{3}*'.'*s^{-1}*'] x'~10^{.(as.character(power))})) + } + } else { + p = p + + ggtitle(bquote(.(type)~' ['*m^{3}*'.'*s^{-1}*'] x'~10^{.(as.character(power))})) + } # if (norm) { @@ -410,8 +209,8 @@ time_panel = function (df_data_code, df_trend_code, type, p_threshold=0.1, missR max(df_data_code$Date)), expand=c(0, 0)) + - scale_y_continuous(breaks=seq(0, maxQ*10, dbrk), - limits=c(0, maxQ*1.1), + scale_y_continuous(breaks=seq(0, maxQN*10, dbrk), + limits=c(0, maxQN*1.1), expand=c(0, 0), labels=label_number(accuracy=accuracy)) @@ -425,7 +224,7 @@ text_panel = function(code, df_meta) { text = paste( "<span style='font-size:18pt'> station <b>", code, "</b></span><br>", "nom : ", df_meta_code$nom, "<br>", - "territoire : ", df_meta_code$territoire, "<br>", + "région hydrographique : ", df_meta_code$region_hydro, "<br>", "position : (", df_meta_code$L93X, "; ", df_meta_code$L93Y, ")", "<br>", "surface : ", df_meta_code$surface_km2, " km<sup>2</sup>", sep='') @@ -440,28 +239,104 @@ text_panel = function(code, df_meta) { -get_color = function (value, min, max, ncolor, palette_name="RdYlBu", reverse=TRUE) { - - if (min == max) { - palette = colorRampPalette(brewer.pal(11, palette_name))(3) - color = palette[2] - return(list(color=color, palette=c(color))) +matrice_panel = function (list_df2plot, df_meta) { + + nbp = length(list_df2plot) + + minTrend = c() + maxTrend = c() + + for (i in 1:nbp) { + + df_trend = list_df2plot[[i]]$trend + p_threshold = list_df2plot[[i]]$p_threshold + + okTrend = df_trend$trend[df_trend$p <= p_threshold] + + minTrend[i] = min(okTrend, na.rm=TRUE) + maxTrend[i] = max(okTrend, na.rm=TRUE) + } + + # Get all different stations code + Code = levels(factor(df_meta$code)) + + Type_mat = c() + Code_mat = c() + Color_mat = c() + + for (code in Code) { + + for (i in 1:nbp) { + df_trend = list_df2plot[[i]]$trend + p_threshold = list_df2plot[[i]]$p_threshold + type = list_df2plot[[i]]$type + + Type_mat[i] = as.character(type) + Code_mat[i] = code + + df_trend_code = df_trend[df_trend$code == code,] + + if (df_trend_code$p <= p_threshold){ + color_res = get_color(df_trend_code$trend, + minTrend[i], + maxTrend[i], + palette_name='perso', + reverse=FALSE) + + color = color_res$color + + } else { + color = 'white' + } + + Color_mat[i] = color + } } + + mat = ggplot() + + geom_tile(aes(x=Type_mat, y=Code_mat, fill=Color_mat)) + + return (mat) +} + + - palette = colorRampPalette(brewer.pal(11, palette_name))(ncolor) + + +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' + ))(ncolor) + + } else { + palette = colorRampPalette(brewer.pal(11, palette_name))(ncolor) + } if (reverse) { palette = rev(palette) } - idNorm = (value - min) / (max - min) - - id = round(idNorm*(ncolor-1) + 1, 0) + palette_cold = palette[1:as.integer(ncolor/2)] + palette_hot = palette[(as.integer(ncolor/2)+1):ncolor] - print(idNorm) - print(id) + ncolor_cold = length(palette_cold) + ncolor_hot = length(palette_hot) - color = palette[id] + if (value < 0) { + idNorm = (value - min) / (0 - min) + id = round(idNorm*(ncolor_cold - 1) + 1, 0) + color = palette_cold[id] + } else { + idNorm = (value - 0) / (max - 0) + id = round(idNorm*(ncolor_hot - 1) + 1, 0) + color = palette_hot[id] + } return(list(color=color, palette=palette)) } @@ -480,3 +355,31 @@ void = ggplot() + geom_blank(aes(1,1)) + axis.ticks = element_blank(), axis.line = element_blank() ) + + + +palette_tester = function () { + + n = 300 + X = 1:n + Y = rep(0, times=n) + + palette = colorRampPalette(c( + '#1a4157', + '#00af9d', + '#fbdd7e', + '#fdb147', + '#fd4659' + ))(n) + + p = ggplot() + + geom_line(aes(x=X, y=Y), color=palette[X], size=10) + + 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() diff --git a/processing/extractBH.R b/processing/extractBH.R index 93e21ef..2ce4f4b 100644 --- a/processing/extractBH.R +++ b/processing/extractBH.R @@ -49,6 +49,40 @@ iQHE = c('0'='qualit '2'='qualité hautes eaux douteuse') +iRegHydro = c('D'='Affluents du Rhin', + 'E'="Fleuves côtiers de l'Artois-Picardie", + 'A'='Rhin', + 'B'='Meuse', + 'F'='Seine aval (Marne incluse)', + 'G'='Fleuves côtiers haut normands', + 'H'='Seine amont', + 'I'='Fleuves côtiers bas normands', + 'J'='Bretagne', + 'K'='Loire', + 'L'='Loire', + 'M'='Loire', + 'N'='Fleuves côtiers au sud de la Loire', + 'O'='Garonne', + 'P'='Dordogne', + 'Q'='Adour', + 'R'='Charente', + 'S'="Fleuves côtiers de l'Adour-Garonne", + 'U'='Saône', + 'V'='Rhône', + 'W'='Isère', + 'X'='Durance', + 'Y'='Fleuves côtiers du Rhône-Méditérannée et Corse', + 'Z'='Îles', + '1'='Guadeloupe', + '2'='Martinique', + '5'='Guyane', + '6'='Guyane', + '7'='Guyane', + '8'='Guyane', + '9'='Guyane', + '4'='Réunion') + + # Get the selection of data from the 'Liste-station_RRSE' file and the BanqueHydro directory get_selection = function (computer_data_path, listdir, listname, cnames=c('code','station', 'BV_km2', 'axe_principal_concerne', 'longueur_serie', 'commentaires', 'choix'), @@ -190,6 +224,8 @@ extractBH_meta = function (computer_data_path, filedir, filename, verbose=TRUE) source='BH' ) + df_meta$region_hydro = iRegHydro[substr(df_meta$code, 1, 1)] + return (df_meta) } else { diff --git a/processing/extractNV.R b/processing/extractNV.R index 8bf6c24..012084e 100644 --- a/processing/extractNV.R +++ b/processing/extractNV.R @@ -3,6 +3,40 @@ library(tools) library(dplyr) +iRegHydro = c('D'='Affluents du Rhin', + 'E'="Fleuves côtiers de l'Artois-Picardie", + 'A'='Rhin', + 'B'='Meuse', + 'F'='Seine aval (Marne incluse)', + 'G'='Fleuves côtiers haut normands', + 'H'='Seine amont', + 'I'='Fleuves côtiers bas normands', + 'J'='Bretagne', + 'K'='Loire', + 'L'='Loire', + 'M'='Loire', + 'N'='Fleuves côtiers au sud de la Loire', + 'O'='Garonne', + 'P'='Dordogne', + 'Q'='Adour', + 'R'='Charente', + 'S'="Fleuves côtiers de l'Adour-Garonne", + 'U'='Saône', + 'V'='Rhône', + 'W'='Isère', + 'X'='Durance', + 'Y'='Fleuves côtiers du Rhône-Méditérannée et Corse', + 'Z'='Îles', + '1'='Guadeloupe', + '2'='Martinique', + '5'='Guyane', + '6'='Guyane', + '7'='Guyane', + '8'='Guyane', + '9'='Guyane', + '4'='Réunion') + + # Extraction of metadata extractNVlist_meta = function (computer_data_path, filedir, listdir, listname, verbose=TRUE) { @@ -57,7 +91,7 @@ extractNVlist_meta = function (computer_data_path, filedir, listdir, listname, v altitude_m=df_meta$Alt, file_path=file.path(dir_path, paste(df_meta$CODE, '.txt', sep='')), - source='NV' + source='NV', ) df_meta = bind_rows(df_meta, @@ -68,6 +102,8 @@ extractNVlist_meta = function (computer_data_path, filedir, listdir, listname, v sep='')))) df_meta = df_meta[order(df_meta$code),] + + df_meta$region_hydro = iRegHydro[substr(df_meta$code, 1, 1)] } else { print(paste('filename', list_path, 'do not exist')) diff --git a/script.R b/script.R index 9950023..93ba19f 100644 --- a/script.R +++ b/script.R @@ -77,7 +77,7 @@ source('processing/extractNV.R', encoding='latin1') source('processing/format.R', encoding='latin1') source('processing/analyse.R', encoding='latin1') source('plotting/panel.R', encoding='latin1') -# source('plotting/layout.R') +source('plotting/layout.R', encoding='latin1') # Usefull library @@ -155,31 +155,31 @@ res_VCN10trend = get_VCN10trend(df_data, df_meta, period) # TIME PANEL # # Plot time panel of debit by stations -# panel(list(df_data, df_data), -# layout_matrix=c(1, 2), -# df_meta=df_meta, -# missRect=list(TRUE, TRUE), -# type=list('Q', 'sqrt(Q)'), -# info_header=TRUE, -# time_header=NULL, -# header_ratio=3, -# figdir=figdir, -# filename_opt='time') - -panel(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, +# panels_layout(list(df_data, df_data), +# layout_matrix=c(1, 2), +# df_meta=df_meta, +# missRect=list(TRUE, TRUE), +# type=list('Q', 'sqrt(Q)'), +# info_header=TRUE, +# time_header=NULL, +# header_ratio=3, +# 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, - header_ratio=2, - figdir=figdir, - filename_opt='') + type=list('Q[A]', 'Q[MNA]', 'V[CN10]'), + missRect=list(TRUE, TRUE, TRUE), + period=period, + info_header=TRUE, + time_header=df_data, + header_ratio=2, + figdir=figdir, + filename_opt='') ### /!\ Removed 185 row(s) containing missing values (geom_path) -> remove NA ### -- GitLab