diff --git a/plotting/datasheet.R b/plotting/datasheet.R index 61ab26389b859679c93080e6eb568de5fc204519..d8ae947504f20c291ed74ff2dd188336b1e55160 100644 --- a/plotting/datasheet.R +++ b/plotting/datasheet.R @@ -147,7 +147,7 @@ datasheet_panel = function (list_df2plot, df_meta, trend_period, info_header, ti df_trend = list_df2plot[[i]]$trend # Extracts the type of the variable type = list_df2plot[[i]]$type - p_threshold = list_df2plot[[i]]$p_threshold + alpha = list_df2plot[[i]]$alpha # Extracts the data corresponding to the code df_data_code = df_data[df_data$code == code,] df_trend_code = df_trend[df_trend$code == code,] @@ -187,7 +187,7 @@ datasheet_panel = function (list_df2plot, df_meta, trend_period, info_header, ti } # If the p value is under the threshold - if (df_trend_code_per$p <= p_threshold) { + if (df_trend_code_per$p <= alpha) { # Stores the mean trend TrendValue_code[j, i, k] = trendValue # Otherwise @@ -338,7 +338,7 @@ datasheet_panel = function (list_df2plot, df_meta, trend_period, info_header, ti # Extracts the trend corresponding to the # current variable df_trend = list_df2plot[[i]]$trend - p_threshold = list_df2plot[[i]]$p_threshold + alpha = list_df2plot[[i]]$alpha unit2day = list_df2plot[[i]]$unit2day missRect = list_df2plot[[i]]$missRect # Extract the variable of the plot @@ -359,7 +359,7 @@ datasheet_panel = function (list_df2plot, df_meta, trend_period, info_header, ti for (j in 1:nPeriod_max) { # If the trend is significant - if (df_trend_code$p[j] <= p_threshold){ + if (df_trend_code$p[j] <= alpha){ # Gets the associated time info Start = tab_Start[k, i, j] End = tab_End[k, i, j] @@ -421,7 +421,7 @@ datasheet_panel = function (list_df2plot, df_meta, trend_period, info_header, ti # Computes the time panel associated to the current variable p = time_panel(df_data_code, df_trend_code, var=var, - type=type, p_threshold=p_threshold, + type=type, alpha=alpha, missRect=missRect, trend_period=trend_period, mean_period=mean_period, axis_xlim=axis_xlim, unit2day=unit2day, grid=FALSE, color=color, @@ -522,7 +522,7 @@ datasheet_panel = function (list_df2plot, df_meta, trend_period, info_header, ti ## 2. OTHER PANEL FOR THE DATASHEET ### 2.1. Time panel -time_panel = function (df_data_code, df_trend_code, var, type, p_threshold=0.1, missRect=FALSE, unit2day=365.25, trend_period=NULL, mean_period=NULL, axis_xlim=NULL, grid=TRUE, ymin_lim=NULL, color=NULL, NspaceMax=NULL, first=FALSE, last=FALSE, lim_pct=10) { +time_panel = function (df_data_code, df_trend_code, var, type, alpha=0.1, missRect=FALSE, unit2day=365.25, trend_period=NULL, mean_period=NULL, axis_xlim=NULL, grid=TRUE, ymin_lim=NULL, color=NULL, NspaceMax=NULL, first=FALSE, last=FALSE, lim_pct=10) { # If 'type' is square root apply it to data if (var == 'sqrt(Q)') { diff --git a/plotting/layout.R b/plotting/layout.R index b382be40746e0d9f92fadec577afa6b35614f6c8..a31cdd89c9e8e25e70a5d7593cebf45585a5b734 100644 --- a/plotting/layout.R +++ b/plotting/layout.R @@ -137,7 +137,7 @@ datasheet_layout = function (df_data, df_meta, layout_matrix, toplot=c('datasheet', 'matrix', 'map'), figdir='', filedir_opt='', filename_opt='', variable='', df_trend=NULL, - p_threshold=0.1, unit2day=365.25, var='', + alpha=0.1, unit2day=365.25, var='', type='', trend_period=NULL, mean_period=NULL, axis_xlim=NULL, missRect=FALSE, time_header=NULL, @@ -192,13 +192,13 @@ datasheet_layout = function (df_data, df_meta, layout_matrix, df_trend = replicate(nbp, df_trend) }} - if (all(class(p_threshold) != 'list')) { - p_threshold = list(p_threshold) + if (all(class(alpha) != 'list')) { + alpha = list(alpha) # If there is only one value - if (length(p_threshold) == 1) { + if (length(alpha) == 1) { # Replicates the value the number of times that there # is of studied variables - p_threshold = replicate(nbp, p_threshold) + alpha = replicate(nbp, alpha) }} # Same @@ -234,7 +234,7 @@ datasheet_layout = function (df_data, df_meta, layout_matrix, # Creates a list that gather all the info for one type of graph df2plot = list(data=df_data[[i]], trend=df_trend[[i]], - p_threshold=p_threshold[[i]], + alpha=alpha[[i]], unit2day=unit2day[[i]], var=var[[i]], type=type[[i]], missRect=missRect[[i]]) @@ -252,14 +252,16 @@ datasheet_layout = function (df_data, df_meta, layout_matrix, # If summarize matrix needs to be plot if ('matrix' %in% toplot) { matrix_panel(list_df2plot, df_meta, trend_period, mean_period, - slice=12, outdirTmp=outdirTmp, A3=TRUE) + slice=19, outdirTmp=outdirTmp, A3=TRUE, + foot_note=foot_note, foot_height=foot_height, resources_path=resources_path, AEAGlogo_file=AEAGlogo_file, INRAElogo_file=INRAElogo_file, FRlogo_file=FRlogo_file,) } # If map needs to be plot if ('map' %in% toplot) { map_panel(list_df2plot, df_meta, - idPer=length(trend_period), + idPer_trend=length(trend_period), + mean_period=mean_period, df_shapefile=df_shapefile, foot_note=foot_note, foot_height=foot_height, diff --git a/plotting/map.R b/plotting/map.R index 174d1d953642c9c4c05b2ba0e1e18868c4974d0f..5c4dcf9b556b4fd5dfa2a88b08f767f6c35e279f 100644 --- a/plotting/map.R +++ b/plotting/map.R @@ -28,7 +28,9 @@ ## 1. MAP PANEL # Generates a map plot of the tendancy of a hydrological variable -map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='', codeLight=NULL, margin=NULL, showSea=TRUE, +map_panel = function (list_df2plot, df_meta, df_shapefile, idPer_trend=1, + mean_period=NULL, outdirTmp='', codeLight=NULL, + margin=NULL, showSea=TRUE, foot_note=FALSE, foot_height=0, resources_path=NULL, AEAGlogo_file=NULL, INRAElogo_file=NULL, @@ -143,7 +145,7 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='' df_trend = list_df2plot[[i]]$trend # Extracts the type of the variable type = list_df2plot[[i]]$type - p_threshold = list_df2plot[[i]]$p_threshold + alpha = list_df2plot[[i]]$alpha # Extracts the data corresponding to the code df_data_code = df_data[df_data$code == code,] # Extracts the trend corresponding to the code @@ -183,7 +185,7 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='' } # If the p value is under the threshold - if (df_trend_code_per$p <= p_threshold){ + if (df_trend_code_per$p <= alpha){ # Stores the mean trend TrendValue_code[j, i, k] = trendValue # Otherwise @@ -198,575 +200,743 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='' minTrendValue = apply(TrendValue_code, c(1, 2), min, na.rm=TRUE) maxTrendValue = apply(TrendValue_code, c(1, 2), max, na.rm=TRUE) + + + + + # If there is a 'mean_period' + if (!is.null(mean_period)) { + # Blank vectors to store info about breaking analysis + Var_mean = c() + Type_mean = c() + Code_mean = c() + DataMean_mean = c() + breakValue_mean = c() + + # Convert 'mean_period' to list + mean_period = as.list(mean_period) + # Number of mean period + nPeriod_mean = length(mean_period) + + # Blank array to store difference of mean between two periods + breakValue_code = array(rep(1, nPeriod_mean*nbp*nCode), + dim=c(nPeriod_mean, nbp, nCode)) + # Blank array to store mean for a temporary period in order + # to compute the difference of mean with a second period + dataMeantmp = array(rep(NA, nbp*nCode), + dim=c(nbp, nCode)) + + # For all period of breaking analysis + for (j in 1:nPeriod_mean) { + # For all the code + for (k in 1:nCode) { + # Gets the code + code = Code[k] + # For all variable + for (i in 1:nbp) { + # Extracts the data corresponding to + # the current variable + df_data = list_df2plot[[i]]$data + # Extract the variable of the plot + var = list_df2plot[[i]]$var + # Extract the type of the variable to plot + type = list_df2plot[[i]]$type + # Extracts the data corresponding to the code + df_data_code = df_data[df_data$code == code,] + + # Get the current start and end of the sub period + Start_mean = mean_period[[j]][1] + End_mean = mean_period[[j]][2] + + # Extract the data corresponding to this sub period + df_data_code_per = + df_data_code[df_data_code$Date >= Start_mean + & df_data_code$Date <= End_mean,] + + # Min max for the sub period + Datemin = min(df_data_code_per$Date) + Datemax = max(df_data_code_per$Date) + + # Mean of the flow over the sub period + dataMean = mean(df_data_code_per$Value, + na.rm=TRUE) + + # If this in not the first period + if (j > 1) { + # Compute the difference of mean + Break = dataMean - dataMeantmp[i, k] + # Otherwise for the first period + } else { + # Stocks NA + Break = NA + } + + # If it is a flow variable + if (type == 'sévérité') { + # Normalises the break by the mean of the + # initial period + breakValue = Break / dataMeantmp[i, k] + # If it is a date variable + } else if (type == 'saisonnalité') { + # Just stocks the break value + breakValue = Break + } + + # Stores the result + breakValue_code[j, i, k] = breakValue + # Stores temporarily the mean of the current period + dataMeantmp[i, k] = dataMean + } + } + } + # Computes the min and the max of the averaged trend for + # all the station + minBreakValue = apply(breakValue_code, c(1, 2), + min, na.rm=TRUE) + maxBreakValue = apply(breakValue_code, c(1, 2), + max, na.rm=TRUE) + } + + if (is.null(mean_period)) { + nPeriod_mean = 1 + } + # Number of ticks for the colorbar nbTick = 10 - # For all variable - for (i in 1:nbp) { - # If there is a specified station code to highlight (mini map) - # and there has already been one loop - if (i > 1 & !is.null(codeLight)) { - # Stop the for loop over the variable - break - } - # Extracts the variable of the plot - var = list_df2plot[[i]]$var - # Extracts the type of variable of the plot - type = list_df2plot[[i]]$type - # Createsa name for the map - outname = paste('map_', var, sep='') - # If there is the verbose option - if (verbose) { - # Prints the name of the map - print(paste('Map for variable : ', var, - " (", round(i/nbp*100, 0), " %)", - sep='')) - } - - # If there is no specified station code to highlight (mini map) - if (is.null(codeLight)) { - # Sets the size of the countour - sizefr = 0.45 - sizebs = 0.4 - sizerv = 0.3 - } else { - sizefr = 0.35 - sizebs = 0.3 - sizerv = 0.2 - } - # Stores the coordonate system - cf = coord_fixed() - # Makes it the default one to remove useless warning - cf$default = TRUE - - # Open a new plot with the personalise theme - map = ggplot() + theme_void() + - # theme(plot.background=element_rect(fill=NA, - # color="#EC4899")) + - # Fixed coordinate system (remove useless warning) - cf + - # Plot the background of France - geom_polygon(data=df_france, - aes(x=long, y=lat, group=group), - color=NA, fill="grey97") - # If the river shapefile exists - if (!is.null(df_river)) { - # Plot the river + for (j in 1:nPeriod_mean) { + # For all variable + for (i in 1:nbp) { + # If there is a specified station code to highlight (mini map) + # and there has already been one loop + if ((i > 1 | j > 1) & !is.null(codeLight)) { + # Stop the for loop over the variable + break + } + # Extracts the variable of the plot + var = list_df2plot[[i]]$var + # Extracts the type of variable of the plot + type = list_df2plot[[i]]$type + + # Createsa name for the map + if (j > 1) { + outname = paste('map_d', var, sep='') + } else { + outname = paste('map_', var, sep='') + } + + n_page = i + nbp*(j-1) + N_page = nbp*nPeriod_mean + # If there is the verbose option + if (verbose) { + # Prints the name of the map + print(paste('Map for variable : ', var, + " (", round(n_page/N_page*100, 0), " %)", + sep='')) + } + + # If there is no specified station code to highlight + # (mini map) + if (is.null(codeLight)) { + # Sets the size of the countour + sizefr = 0.45 + sizebs = 0.4 + sizerv = 0.3 + } else { + sizefr = 0.35 + sizebs = 0.3 + sizerv = 0.2 + } + + # Stores the coordonate system + cf = coord_fixed() + # Makes it the default one to remove useless warning + cf$default = TRUE + + # Open a new plot with the personalise theme + map = ggplot() + theme_void() + + # theme(plot.background=element_rect(fill=NA, + # color="#EC4899")) + + # Fixed coordinate system (remove useless warning) + cf + + # Plot the background of France + geom_polygon(data=df_france, + aes(x=long, y=lat, group=group), + color=NA, fill="grey97") + # If the river shapefile exists + if (!is.null(df_river)) { + # Plot the river + map = map + + geom_path(data=df_river, + aes(x=long, y=lat, group=group), + color="grey85", size=sizerv) + } + map = map + - geom_path(data=df_river, - aes(x=long, y=lat, group=group), - color="grey85", size=sizerv) - } - - map = map + - # Plot the hydrological basin - geom_polygon(data=df_bassin, - aes(x=long, y=lat, group=group), - color="grey70", fill=NA, size=sizebs) + - # Plot the hydrological sub-basin - geom_polygon(data=df_subbassin, - aes(x=long, y=lat, group=group), - color="grey70", fill=NA, size=sizebs) + - # Plot the countour of France - geom_polygon(data=df_france, - aes(x=long, y=lat, group=group), - color="grey40", fill=NA, size=sizefr) - - # If the sea needs to be shown - if (showSea) { - # Leaves space around the France - xlim = c(295000, 790000) - ylim = c(6125000, 6600000) - # Otherwise - } else { - # Leaves minimal space around France - xlim = c(305000, 790000) - ylim = c(6135000, 6600000) - } + # Plot the hydrological basin + geom_polygon(data=df_bassin, + aes(x=long, y=lat, group=group), + color="grey70", fill=NA, size=sizebs) + + # Plot the hydrological sub-basin + geom_polygon(data=df_subbassin, + aes(x=long, y=lat, group=group), + color="grey70", fill=NA, size=sizebs) + + # Plot the countour of France + geom_polygon(data=df_france, + aes(x=long, y=lat, group=group), + color="grey40", fill=NA, size=sizefr) + + # If the sea needs to be shown + if (showSea) { + # Leaves space around the France + xlim = c(295000, 790000) + ylim = c(6125000, 6600000) + # Otherwise + } else { + # Leaves minimal space around France + xlim = c(305000, 790000) + ylim = c(6135000, 6600000) + } - # If there is no specified station code to highlight (mini map) - if (is.null(codeLight)) { - # Sets a legend scale start - xmin = gpct(4, xlim, shift=TRUE) - # Sets graduations - xint = c(0, 10*1E3, 50*1E3, 100*1E3) - # Sets the y postion - ymin = gpct(5, ylim, shift=TRUE) - # Sets the height of graduations - ymax = ymin + gpct(1, ylim) - # Size of the value - size = 3 - # Size of the 'km' unit - sizekm = 2.5 - # If there is a specified station code - } else { - # Same but with less graduation and smaller size - xmin = gpct(2, xlim, shift=TRUE) - xint = c(0, 100*1E3) - ymin = gpct(1, ylim, shift=TRUE) - ymax = ymin + gpct(3, ylim) - size = 2 - sizekm = 1.5 - } - - map = map + - # Adds the base line of the scale - geom_line(aes(x=c(xmin, max(xint)+xmin), - y=c(ymin, ymin)), - color="grey40", size=0.2) + - # Adds the 'km' unit - annotate("text", - x=max(xint)+xmin+gpct(1, xlim), y=ymin, - vjust=0, hjust=0, label="km", - color="grey40", size=sizekm) - # For all graduations - for (x in xint) { + # If there is no specified station code to highlight (mini map) + if (is.null(codeLight)) { + # Sets a legend scale start + xmin = gpct(4, xlim, shift=TRUE) + # Sets graduations + xint = c(0, 10*1E3, 50*1E3, 100*1E3) + # Sets the y postion + ymin = gpct(5, ylim, shift=TRUE) + # Sets the height of graduations + ymax = ymin + gpct(1, ylim) + # Size of the value + size = 3 + # Size of the 'km' unit + sizekm = 2.5 + # If there is a specified station code + } else { + # Same but with less graduation and smaller size + xmin = gpct(2, xlim, shift=TRUE) + xint = c(0, 100*1E3) + ymin = gpct(1, ylim, shift=TRUE) + ymax = ymin + gpct(3, ylim) + size = 2 + sizekm = 1.5 + } + map = map + - # Draws the tick - annotate("segment", - x=x+xmin, xend=x+xmin, y=ymin, yend=ymax, - color="grey40", size=0.2) + - # Adds the value + # Adds the base line of the scale + geom_line(aes(x=c(xmin, max(xint)+xmin), + y=c(ymin, ymin)), + color="grey40", size=0.2) + + # Adds the 'km' unit annotate("text", - x=x+xmin, y=ymax+gpct(0.5, ylim), - vjust=0, hjust=0.5, label=x/1E3, - color="grey40", size=size) - } - - map = map + - # Allows to crop shapefile without graphical problem - coord_sf(xlim=xlim, ylim=ylim, - expand=FALSE) - - # If there is no margins specified - if (is.null(margin)) { - # Sets all margins to 0 - map = map + - theme(plot.margin=margin(t=0, r=0, b=0, l=0, unit="mm")) - # Otherwise - } else { - # Sets margins to the given ones - map = map + - theme(plot.margin=margin) - } + x=max(xint)+xmin+gpct(1, xlim), y=ymin, + vjust=0, hjust=0, label="km", + color="grey40", size=sizekm) + # For all graduations + for (x in xint) { + map = map + + # Draws the tick + annotate("segment", + x=x+xmin, xend=x+xmin, y=ymin, yend=ymax, + color="grey40", size=0.2) + + # Adds the value + annotate("text", + x=x+xmin, y=ymax+gpct(0.5, ylim), + vjust=0, hjust=0.5, label=x/1E3, + color="grey40", size=size) + } + + map = map + + # Allows to crop shapefile without graphical problem + coord_sf(xlim=xlim, ylim=ylim, + expand=FALSE) + + # If there is no margins specified + if (is.null(margin)) { + # Sets all margins to 0 + map = map + + theme(plot.margin=margin(t=0, r=0, b=0, l=0, + unit="mm")) + # Otherwise + } else { + # Sets margins to the given ones + map = map + + theme(plot.margin=margin) + } - # Blank vector to store data about station - lon = c() - lat = c() - fill = c() - shape = c() - trend = c() - p_threshold_Ok = c() - # For all code - for (k in 1:nCode) { - # Gets the code - code = Code[k] - # Extracts the data corresponding to the current variable - df_data = list_df2plot[[i]]$data - # Extracts the trend corresponding to the - # current variable - df_trend = list_df2plot[[i]]$trend - p_threshold = list_df2plot[[i]]$p_threshold - # Extracts the data corresponding to the code - df_data_code = df_data[df_data$code == code,] - # Extracts the trend corresponding to the code - df_trend_code = df_trend[df_trend$code == code,] + # Blank vector to store data about station + lon = c() + lat = c() + fill = c() + shape = c() + Value = c() + alpha_Ok = c() + # For all code + for (k in 1:nCode) { + # Gets the code + code = Code[k] + + if (j > 1) { + value = breakValue_code[j, i, k] + minValue = minBreakValue[j, i] + maxValue = maxBreakValue[j, i] + pvalue = 0 + + } else { + # Extracts the data corresponding to the + # current variable + df_data = list_df2plot[[i]]$data + # Extracts the trend corresponding to the + # current variable + df_trend = list_df2plot[[i]]$trend + # Gets the risk of the test + alpha = list_df2plot[[i]]$alpha + # Extracts the data corresponding to the code + df_data_code = df_data[df_data$code == code,] + # Extracts the trend corresponding to the code + df_trend_code = df_trend[df_trend$code == code,] + + # Gets the associated time info + Start = tab_Start[k, i, idPer_trend] + End = tab_End[k, i, idPer_trend] + Periods = tab_Periods[k, i, idPer_trend] + + # Extracts the corresponding data for the period + df_data_code_per = + df_data_code[df_data_code$Date >= Start + & df_data_code$Date <= End,] + # Same for trend + df_trend_code_per = + df_trend_code[df_trend_code$period_start == Start + & df_trend_code$period_end == End,] + + # Computes the number of trend analysis selected + Ntrend = nrow(df_trend_code_per) + # If there is more than one trend on the same period + if (Ntrend > 1) { + # Takes only the first because they are similar + df_trend_code_per = df_trend_code_per[1,] + } + + # If it is a flow variable + if (type == 'sévérité') { + # Computes the mean of the data on the period + dataMean = mean(df_data_code_per$Value, + na.rm=TRUE) + # Normalises the trend value by the mean + # of the data + value = df_trend_code_per$trend / dataMean + # If it is a date variable + } else if (type == 'saisonnalité') { + value = df_trend_code_per$trend + } + + minValue = minTrendValue[idPer_trend, i] + maxValue = maxTrendValue[idPer_trend, i] + pvalue = df_trend_code_per$p + + } + + # Computes the color associated to the mean trend + color_res = get_color(value, + minValue, + maxValue, + palette_name='perso', + reverse=TRUE, + ncolor=256) + # Computes the colorbar info + palette_res = get_palette(minValue, + maxValue, + palette_name='perso', + reverse=TRUE, + ncolor=256, + nbTick=nbTick) + + # If it is significative + if (pvalue <= alpha){ + # The computed color is stored + filltmp = color_res + # If the mean tend is positive + if (value >= 0) { + # Uses a triangle up for the shape + # of the marker + shapetmp = 24 + # If negative + } else { + # Uses a triangle down for the shape + # of the marker + shapetmp = 25 + } + # If it is not significative + } else { + # The fill color is grey + filltmp = 'grey97' + # The marker is a circle + shapetmp = 21 + } - # Gets the associated time info - # Start = Start_code[Code_code == code][[1]][idPer] - # End = End_code[Code_code == code][[1]][idPer] - - # Gets the associated time info - Start = tab_Start[k, i, idPer] - End = tab_End[k, i, idPer] - Periods = tab_Periods[k, i, idPer] - - # Extracts the corresponding data for the period - df_data_code_per = - df_data_code[df_data_code$Date >= Start - & df_data_code$Date <= End,] - # Same for trend - df_trend_code_per = - df_trend_code[df_trend_code$period_start == Start - & df_trend_code$period_end == End,] - - # Computes the number of trend analysis selected - Ntrend = nrow(df_trend_code_per) - # If there is more than one trend on the same period - if (Ntrend > 1) { - # Takes only the first because they are similar - df_trend_code_per = df_trend_code_per[1,] + # Extracts the localisation of the current station + lontmp = + df_meta[df_meta$code == code,]$L93X_m_BH + lattmp = + df_meta[df_meta$code == code,]$L93Y_m_BH + + # Stores all the parameters + lon = c(lon, lontmp) + lat = c(lat, lattmp) + fill = c(fill, filltmp) + shape = c(shape, shapetmp) + Value = c(Value, value) + # If the trend analysis is significative a TRUE is stored + alpha_Ok = c(alpha_Ok, + pvalue <= alpha) } + # Creates a tibble to stores all the data to plot + plot_map = tibble(lon=lon, lat=lat, fill=fill, + shape=shape, code=Code) + + # If there is no specified station code to highlight + # (mini map) + if (is.null(codeLight)) { + map = map + + # Plots the trend point + geom_point(data=plot_map, + aes(x=lon, y=lat), + shape=shape, size=5, stroke=1, + color='grey50', fill=fill) + # If there is a specified station code + } else { + # Extract data of all stations not to highlight + plot_map_codeNo = plot_map[plot_map$code != codeLight,] + # Extract data of the station to highlight + plot_map_code = plot_map[plot_map$code == codeLight,] + + # Plots only the localisation + map = map + + # For all stations not to highlight + geom_point(data=plot_map_codeNo, + aes(x=lon, y=lat), + shape=21, size=0.5, stroke=0.5, + color='grey50', fill='grey50') + + # For the station to highlight + geom_point(data=plot_map_code, + aes(x=lon, y=lat), + shape=21, size=1.5, stroke=0.5, + color='#00A3A8', fill='#00A3A8') + } + + # Extracts the position of the tick of the colorbar + posTick = palette_res$posTick + # Extracts the label of the tick of the colorbar + labTick = palette_res$labTick + # Extracts the color corresponding to the tick of the colorbar + colTick = palette_res$colTick + + # Spreading of the colorbar + valNorm = nbTick * 10 + # Normalisation of the position of ticks + ytick = posTick / max(posTick) * valNorm # If it is a flow variable if (type == 'sévérité') { - # Computes the mean of the data on the period - dataMean = mean(df_data_code_per$Value, na.rm=TRUE) - # Normalises the trend value by the mean of the data - trendValue = df_trend_code_per$trend / dataMean - # If it is a date variable + # Formatting of label in pourcent + labTick = as.character(round(labTick*100, 2)) + # If it is a date variable } else if (type == 'saisonnalité') { - trendValue = df_trend_code_per$trend + # Formatting of label + labTick = as.character(round(labTick, 2)) } - - # Computes the color associated to the mean trend - color_res = get_color(trendValue, - minTrendValue[idPer, i], - maxTrendValue[idPer, i], - palette_name='perso', - reverse=TRUE, - ncolor=256) - # Computes the colorbar info - palette_res = get_palette(minTrendValue[idPer, i], - maxTrendValue[idPer, i], - palette_name='perso', - reverse=TRUE, - ncolor=256, - nbTick=nbTick) - - # If it is significative - if (df_trend_code_per$p <= p_threshold){ - # The computed color is stored - filltmp = color_res - # If the mean tend is positive - if (trendValue >= 0) { - # Uses a triangle up for the shape of the marker - shapetmp = 24 - # If negative - } else { - # Uses a triangle down for the shape of the marker - shapetmp = 25 + + # X position of ticks all similar + xtick = rep(0, times=nbTick) + + # Creates a tibble to store all parameters of colorbar + plot_palette = tibble(xtick=xtick, ytick=ytick, + colTick=colTick, labTick=labTick) + + # New plot with void theme + title = ggplot() + theme_void() + + # Plots separation line + geom_line(aes(x=c(-0.3, 3.7), y=c(0.05, 0.05)), + size=0.6, color="#00A3A8") + + # Writes title + geom_shadowtext(data=tibble(x=-0.3, y=0.2, + label=var), + aes(x=x, y=y, label=label), + fontface="bold", + color="#00A3A8", + bg.colour="white", + hjust=0, vjust=0, size=10) + + # X axis + scale_x_continuous(limits=c(-0.3, 1 + 3), + expand=c(0, 0)) + + # Y axis + scale_y_continuous(limits=c(0, 10), + expand=c(0, 0)) + + # Margin + theme(plot.margin=margin(t=0, r=0, b=0, l=0, unit="mm")) + + # New plot with void theme + pal = ggplot() + theme_void() + + # Plots the point of the colorbar + geom_point(data=plot_palette, + aes(x=xtick, y=ytick), + shape=21, size=5, stroke=1, + color='white', fill=colTick) + + if (j > 1) { + ValueName = "Écarts observées" + # If it is a flow variable + if (type == 'sévérité') { + unit = bquote(bold("(%)")) + # If it is a date variable + } else if (type == 'saisonnalité') { + unit = bquote(bold("(jour)")) } - # If it is not significative } else { - # The fill color is grey - filltmp = 'grey97' - # The marker is a circle - shapetmp = 21 + ValueName = "Tendances observées" + # If it is a flow variable + if (type == 'sévérité') { + unit = bquote(bold("(% par an)")) + # If it is a date variable + } else if (type == 'saisonnalité') { + unit = bquote(bold("(jour par an)")) + } + } + + pal = pal + + # Name of the colorbar + annotate('text', + x=-0.3, y= valNorm + 23, + label=ValueName, + hjust=0, vjust=0.5, + size=6, color='grey40') + + # Unit legend of the colorbar + annotate('text', + x=-0.2, y= valNorm + 13, + label=unit, + hjust=0, vjust=0.5, + size=4, color='grey40') + # For all the ticks + for (id in 1:nbTick) { + pal = pal + + # Adds the value + annotate('text', x=xtick[id]+0.3, + y=ytick[id], + label=bquote(bold(.(labTick[id]))), + hjust=0, vjust=0.7, + size=3, color='grey40') } - # Extracts the localisation of the current station - lontmp = - df_meta[df_meta$code == code,]$L93X_m_BH - lattmp = - df_meta[df_meta$code == code,]$L93Y_m_BH - - # Stores all the parameters - lon = c(lon, lontmp) - lat = c(lat, lattmp) - fill = c(fill, filltmp) - shape = c(shape, shapetmp) - trend = c(trend, trendValue) - # If the trend analysis is significative a TRUE is stored - p_threshold_Ok = c(p_threshold_Ok, - df_trend_code_per$p <= p_threshold) - } - # Creates a tibble to stores all the data to plot - plot_map = tibble(lon=lon, lat=lat, fill=fill, - shape=shape, code=Code) + yUp = -20 + yNone = -29 - # If there is no specified station code to highlight (mini map) - if (is.null(codeLight)) { - map = map + - # Plots the trend point - geom_point(data=plot_map, - aes(x=lon, y=lat), - shape=shape, size=5, stroke=1, - color='grey50', fill=fill) - # If there is a specified station code - } else { - # Extract data of all stations not to highlight - plot_map_codeNo = plot_map[plot_map$code != codeLight,] - # Extract data of the station to highlight - plot_map_code = plot_map[plot_map$code == codeLight,] - - # Plots only the localisation - map = map + - # For all stations not to highlight - geom_point(data=plot_map_codeNo, - aes(x=lon, y=lat), - shape=21, size=0.5, stroke=0.5, - color='grey50', fill='grey50') + - # For the station to highlight - geom_point(data=plot_map_code, - aes(x=lon, y=lat), - shape=21, size=1.5, stroke=0.5, - color='#00A3A8', fill='#00A3A8') - } - - # Extracts the position of the tick of the colorbar - posTick = palette_res$posTick - # Extracts the label of the tick of the colorbar - labTick = palette_res$labTick - # Extracts the color corresponding to the tick of the colorbar - colTick = palette_res$colTick - - # Spreading of the colorbar - valNorm = nbTick * 10 - # Normalisation of the position of ticks - ytick = posTick / max(posTick) * valNorm - - # If it is a flow variable - if (type == 'sévérité') { - # Formatting of label in pourcent - labTick = as.character(round(labTick*100, 2)) - # If it is a date variable - } else if (type == 'saisonnalité') { - # Formatting of label - labTick = as.character(round(labTick, 2)) - } - - # X position of ticks all similar - xtick = rep(0, times=nbTick) - - # Creates a tibble to store all parameters of colorbar - plot_palette = tibble(xtick=xtick, ytick=ytick, - colTick=colTick, labTick=labTick) - - # New plot with void theme - title = ggplot() + theme_void() + - # Plots separation line - geom_line(aes(x=c(-0.3, 3.7), y=c(0.05, 0.05)), - size=0.6, color="#00A3A8") + - # Writes title - geom_shadowtext(data=tibble(x=-0.3, y=0.2, - label=var), - aes(x=x, y=y, label=label), - fontface="bold", - color="#00A3A8", - bg.colour="white", - hjust=0, vjust=0, size=10) + - # X axis - scale_x_continuous(limits=c(-0.3, 1 + 3), - expand=c(0, 0)) + - # Y axis - scale_y_continuous(limits=c(0, 10), - expand=c(0, 0)) + - # Margin - theme(plot.margin=margin(t=0, r=0, b=0, l=0, unit="mm")) - - # New plot with void theme - pal = ggplot() + theme_void() + - # Plots the point of the colorbar - geom_point(data=plot_palette, - aes(x=xtick, y=ytick), - shape=21, size=5, stroke=1, - color='white', fill=colTick) - - # If it is a flow variable - if (type == 'sévérité') { - unit = bquote(bold("% par an")) - # If it is a date variable - } else if (type == 'saisonnalité') { - unit = bquote(bold("jour par an")) - } - - pal = pal + - # Name of the colorbar - annotate('text', - x=-0.3, y= valNorm + 23, - label="Tendance observée", - hjust=0, vjust=0.5, - size=6, color='grey40') + - # Unit legend of the colorbar - annotate('text', - x=-0.2, y= valNorm + 13, - label=unit, - hjust=0, vjust=0.5, - size=4, color='grey40') - # For all the ticks - for (j in 1:nbTick) { + if (j > 1) { + upLabel = bquote(bold("Hausse")) + noneLabel = NULL + downLabel = bquote(bold("Baisse")) + yDown = -29 + } else { + upLabel = bquote(bold("Hausse significative à 10%")) + noneLabel = bquote(bold("Non significatif à 10%")) + downLabel = bquote(bold("Baisse significative à 10%")) + yDown = -40 + } + pal = pal + - # Adds the value - annotate('text', x=xtick[j]+0.3, - y=ytick[j], - label=bquote(bold(.(labTick[j]))), - hjust=0, vjust=0.7, + # Up triangle in the marker legend + geom_point(aes(x=0, y=yUp), + shape=24, size=4, stroke=1, + color='grey50', fill='grey97') + + # Up triangle text legend + annotate('text', + x=0.3, y=yUp, + label=upLabel, + hjust=0, vjust=0.5, size=3, color='grey40') - } - pal = pal + - # Up triangle in the marker legend - geom_point(aes(x=0, y=-20), - shape=24, size=4, stroke=1, - color='grey50', fill='grey97') + - # Up triangle text legend - annotate('text', - x=0.3, y=-20, - label=bquote(bold("Hausse significative à 10%")), - hjust=0, vjust=0.5, - size=3, color='grey40') - - pal = pal + - # Circle in the marker legend - geom_point(aes(x=0, y=-29), - shape=21, size=4, stroke=1, - color='grey50', fill='grey97') + - # Circle text legend - annotate('text', - x=0.3, y=-29, - label=bquote(bold("Non significatif à 10%")), - hjust=0, vjust=0.7, - size=3, color='grey40') - - pal = pal + - # Down triangle in the marker legend - geom_point(aes(x=0, y=-40), - shape=25, size=4, stroke=1, - color='grey50', fill='grey97') + - # Down triangle text legend - annotate('text', - x=0.3, y=-40, - label=bquote(bold("Baisse significative à 10%")), - hjust=0, vjust=0.5, - size=3, color='grey40') - - # Normalises all the trend values for each station - # according to the colorbar - yTrend = (trend - minTrendValue[idPer, i]) / - (maxTrendValue[idPer, i] - minTrendValue[idPer, i]) * valNorm - # Takes only the significative ones - yTrend = yTrend[p_threshold_Ok] - - ## Random distribution ## - # xTrend = rnorm(length(yTrend), mean=1.75, sd=0.1) - - ## Histogram distribution ## - # Computes the histogram of the trend - res_hist = hist(yTrend, breaks=ytick, plot=FALSE) - # Extracts the number of counts per cells - counts = res_hist$counts - # Extracts limits of cells - breaks = res_hist$breaks - # Extracts middle of cells - mids = res_hist$mids - - # Blank vectors to store position of points of - # the distribution to plot - xTrend = c() - yTrend = c() - # Start X position of the distribution - start_hist = 1.25 - # X separation bewteen point - hist_sep = 0.15 - # For all cells of the histogram - for (ii in 1:length(mids)) { - # If the count in the current cell is not zero - if (counts[ii] != 0) { - # Stores the X positions of points of the distribution - # for the current cell - xTrend = c(xTrend, - seq(start_hist, - start_hist+(counts[ii]-1)*hist_sep, - by=hist_sep)) + if (!is.null(noneLabel)) { + pal = pal + + # Circle in the marker legend + geom_point(aes(x=0, y=yNone), + shape=21, size=4, stroke=1, + color='grey50', fill='grey97') + + # Circle text legend + annotate('text', + x=0.3, y=yNone, + label=noneLabel, + hjust=0, vjust=0.7, + size=3, color='grey40') } - # Stores the Y position which is the middle of the - # current cell the number of times it has been counted - yTrend = c(yTrend, rep(mids[ii], times=counts[ii])) - } - - # Makes a tibble to plot the trend distribution - plot_trend = tibble(xTrend=xTrend, yTrend=yTrend) - - pal = pal + - # Plots the point of the trend distribution - geom_point(data=plot_trend, - aes(x=xTrend, y=yTrend), - # shape=21, size=1, - # color="grey20", fill="grey20") - alpha=0.4) - - if (type == 'sévérité') { - labelArrow = 'Plus sévère' - } else if (type == 'saisonnalité') { - labelArrow = 'Plus tôt' - } + + pal = pal + + # Down triangle in the marker legend + geom_point(aes(x=0, y=yDown), + shape=25, size=4, stroke=1, + color='grey50', fill='grey97') + + # Down triangle text legend + annotate('text', + x=0.3, y=yDown, + label=downLabel, + hjust=0, vjust=0.5, + size=3, color='grey40') + + # Normalises all the trend values for each station + # according to the colorbar + if (j > 1) { + yValue = (Value - minBreakValue[j, i]) / (maxBreakValue[j, i] - minBreakValue[j, i]) * valNorm + } else { + yValue = (Value - minTrendValue[idPer_trend, i]) / (maxTrendValue[idPer_trend, i] - minTrendValue[idPer_trend, i]) * valNorm + } + + # Takes only the significative ones + yValue = yValue[alpha_Ok] + + # Histogram distribution + # Computes the histogram of the trend + res_hist = hist(yValue, breaks=ytick, plot=FALSE) + # Extracts the number of counts per cells + counts = res_hist$counts + # Extracts limits of cells + breaks = res_hist$breaks + # Extracts middle of cells + mids = res_hist$mids + + # Blank vectors to store position of points of + # the distribution to plot + xValue = c() + yValue = c() + # Start X position of the distribution + start_hist = 1.25 + # X separation bewteen point + hist_sep = 0.15 + # For all cells of the histogram + for (ii in 1:length(mids)) { + # If the count in the current cell is not zero + if (counts[ii] != 0) { + # Stores the X positions of points of the distribution + # for the current cell + xValue = c(xValue, + seq(start_hist, + start_hist+(counts[ii]-1)*hist_sep, + by=hist_sep)) + } + # Stores the Y position which is the middle of the + # current cell the number of times it has been counted + yValue = c(yValue, rep(mids[ii], times=counts[ii])) + } + + # Makes a tibble to plot the trend distribution + plot_value = tibble(xValue=xValue, yValue=yValue) + + pal = pal + + # Plots the point of the trend distribution + geom_point(data=plot_value, + aes(x=xValue, y=yValue), + # shape=21, size=1, + # color="grey20", fill="grey20") + alpha=0.4) - xArrow = 3.2 - - pal = pal + - # Arrow to show a worsening of the situation - geom_segment(aes(x=xArrow, y=valNorm*0.75, - xend=xArrow, yend=valNorm*0.25), - color='grey50', size=0.3, - arrow=arrow(length=unit(2, "mm"))) + - # Text associated to the arrow - annotate('text', - x=xArrow+0.1, y=valNorm*0.5, - label=labelArrow, - angle=90, - hjust=0.5, vjust=1, - size=3, color='grey50') - - pal = pal + - # X axis of the colorbar - scale_x_continuous(limits=c(-0.3, 1 + 3), - expand=c(0, 0)) + - # Y axis of the colorbar - scale_y_continuous(limits=c(-60, valNorm + 35), - expand=c(0, 0)) + - # Margin of the colorbar - theme(plot.margin=margin(t=0, r=0, b=0, l=0, unit="mm")) - - # If there is a foot note - if (foot_note) { - foot = foot_panel('carte de tendance observée', - i, nbp, resources_path, - AEAGlogo_file, INRAElogo_file, - FRlogo_file, foot_height) - - # Stores the map, the title and the colorbar in a list - P = list(map, title, pal, foot) - LM = matrix(c(1, 1, 1, 2, - 1, 1, 1, 3, - 4, 4, 4, 4), - nrow=3, byrow=TRUE) - } else { - foot_height = 0 - # Stores the map, the title and the colorbar in a list - P = list(map, title, pal) - LM = matrix(c(1, 1, 1, 2, - 1, 1, 1, 3), - nrow=2, byrow=TRUE) - } - - LMcol = ncol(LM) - LMrow = nrow(LM) - - LM = rbind(rep(99, times=LMcol), LM, rep(99, times=LMcol)) - LMrow = nrow(LM) - LM = cbind(rep(99, times=LMrow), LM, rep(99, times=LMrow)) - LMcol = ncol(LM) - - margin_height = 0.5 - height = 21 - width = 29.7 + if (type == 'sévérité') { + labelArrow = 'Plus sévère' + } else if (type == 'saisonnalité') { + labelArrow = 'Plus tôt' + } - row_height = (height - 2*margin_height - foot_height) / (LMrow - 3) + xArrow = 3.2 + + pal = pal + + # Arrow to show a worsening of the situation + geom_segment(aes(x=xArrow, y=valNorm*0.75, + xend=xArrow, yend=valNorm*0.25), + color='grey50', size=0.3, + arrow=arrow(length=unit(2, "mm"))) + + # Text associated to the arrow + annotate('text', + x=xArrow+0.1, y=valNorm*0.5, + label=labelArrow, + angle=90, + hjust=0.5, vjust=1, + size=3, color='grey50') + + pal = pal + + # X axis of the colorbar + scale_x_continuous(limits=c(-0.3, 1 + 3), + expand=c(0, 0)) + + # Y axis of the colorbar + scale_y_continuous(limits=c(-60, valNorm + 35), + expand=c(0, 0)) + + # Margin of the colorbar + theme(plot.margin=margin(t=0, r=0, b=0, l=0, unit="mm")) - Hcut = LM[, 2] - heightLM = rep(row_height, times=LMrow) - heightLM[Hcut == 4] = foot_height - heightLM[Hcut == 99] = margin_height + if (j > 1) { + footName = 'carte des écarts observés' + } else { + footName = 'carte des tendances observées' + } + + # If there is a foot note + if (foot_note) { + foot = foot_panel(footName, + n_page, N_page, resources_path, + AEAGlogo_file, INRAElogo_file, + FRlogo_file, foot_height) + + # Stores the map, the title and the colorbar in a list + P = list(map, title, pal, foot) + LM = matrix(c(1, 1, 1, 2, + 1, 1, 1, 3, + 4, 4, 4, 4), + nrow=3, byrow=TRUE) + } else { + foot_height = 0 + # Stores the map, the title and the colorbar in a list + P = list(map, title, pal) + LM = matrix(c(1, 1, 1, 2, + 1, 1, 1, 3), + nrow=2, byrow=TRUE) + } + id_foot = 4 + + LMcol = ncol(LM) + LMrow = nrow(LM) + + LM = rbind(rep(99, times=LMcol), LM, rep(99, times=LMcol)) + LMrow = nrow(LM) + LM = cbind(rep(99, times=LMrow), LM, rep(99, times=LMrow)) + LMcol = ncol(LM) + + margin_height = 0.5 + height = 21 + width = 29.7 - col_width = (width - 2*margin_height) / (LMcol - 2) - - Wcut = LM[(nrow(LM)-1),] - widthLM = rep(col_width, times=LMcol) - widthLM[Wcut == 99] = margin_height - - # Arranges the graphical object - plot = grid.arrange(grobs=P, layout_matrix=LM, - heights=heightLM, widths=widthLM) - - - # If there is no specified station code to highlight (mini map) - if (is.null(codeLight)) { - # Saving matrix plot - ggsave(plot=plot, - path=outdirTmp, - filename=paste(outname, '.pdf', sep=''), - width=width, height=height, units='cm', dpi=100) + row_height = (height - 2*margin_height - foot_height) / (LMrow - 3) + + Hcut = LM[, 2] + heightLM = rep(row_height, times=LMrow) + heightLM[Hcut == id_foot] = foot_height + heightLM[Hcut == 99] = margin_height + + col_width = (width - 2*margin_height) / (LMcol - 2) + + Wcut = LM[(nrow(LM)-1),] + widthLM = rep(col_width, times=LMcol) + widthLM[Wcut == 99] = margin_height + + # Arranges the graphical object + plot = grid.arrange(grobs=P, layout_matrix=LM, + heights=heightLM, widths=widthLM) + + + # If there is no specified station code to highlight (mini map) + if (is.null(codeLight)) { + # Saving matrix plot + ggsave(plot=plot, + path=outdirTmp, + filename=paste(outname, '.pdf', sep=''), + width=width, height=height, units='cm', dpi=100) + } } } # Returns the map object diff --git a/plotting/matrix.R b/plotting/matrix.R index 84bebd10dbce3fef1ffe062fb5bc2ca3890cb668..679482fd716badb97e441763ba6a5b43f6302403 100644 --- a/plotting/matrix.R +++ b/plotting/matrix.R @@ -28,7 +28,11 @@ ## 1. MATRIX PANEL # Generates a summarizing matrix of the trend analyses of all station for different hydrological variables and periods. Also shows difference of means between specific periods. -matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice=NULL, outdirTmp='', outnameTmp='matrix', title=NULL, A3=FALSE) { +matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice=NULL, outdirTmp='', outnameTmp='matrix', title=NULL, A3=FALSE, + foot_note=FALSE, + foot_height=0, resources_path=NULL, + AEAGlogo_file=NULL, INRAElogo_file=NULL, + FRlogo_file=NULL) { # Number of variable/plot nbp = length(list_df2plot) @@ -140,7 +144,7 @@ matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice df_trend = list_df2plot[[i]]$trend # Extracts the type of the variable type = list_df2plot[[i]]$type - p_threshold = list_df2plot[[i]]$p_threshold + alpha = list_df2plot[[i]]$alpha # Extracts the data corresponding to the code df_data_code = df_data[df_data$code == code,] # Extracts the trend corresponding to the code @@ -181,7 +185,7 @@ matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice } # If the p value is under the threshold - if (df_trend_code_per$p <= p_threshold){ + if (df_trend_code_per$p <= alpha){ # Stores the averaged trend TrendValue_code[j, i, k] = trendValue # Otherwise @@ -221,7 +225,7 @@ matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice # Extracts the trend corresponding to the # current variable df_trend = list_df2plot[[i]]$trend - p_threshold = list_df2plot[[i]]$p_threshold + alpha = list_df2plot[[i]]$alpha # Extract the variable of the plot var = list_df2plot[[i]]$var # Extract the type of the variable to plot @@ -267,7 +271,7 @@ matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice } # If the p value is under the threshold - if (df_trend_code_per$p <= p_threshold){ + if (df_trend_code_per$p <= alpha){ # Gets the color associated to the averaged trend color_res = get_color(trendValue, minTrendValue[j, i], @@ -406,7 +410,8 @@ matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice minBreakValue = apply(breakValue_code, c(1, 2), min, na.rm=TRUE) maxBreakValue = apply(breakValue_code, c(1, 2), - max, na.rm=TRUE) + max, na.rm=TRUE) + # Blanks vector to store color info Fill_mean = c() Color_mean = c() @@ -446,8 +451,17 @@ matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice slice = nCode } + allType = c() + for (i in 1:nbp) { + allType = c(allType, list_df2plot[[i]]$type) + } + + countType = rle(sort(allType)) + df_countType = tibble(type=countType$values, n=countType$lengths) + nbpMax = max(df_countType$n) + # Gets all the different type of plots - Type = levels(factor(Type_trend)) + Type = levels(factor(allType)) nbType = length(Type) # For all the type of plots for (itype in 1:nbType) { @@ -471,12 +485,15 @@ matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice nMat = as.integer(nsubCodefL/slice) + 1 # For all the pages for (iMat in 1:nMat) { + n_page = ifL + nfL*(itype-1) + N_page = nfL*2 + # Print the matrix name print(paste('Matrix ', iMat, '/', nMat, ' of ', type, ' for region : ', fL, " (", - round((ifL + nfL*(itype-1)) / (nfL*2) * 100, + round(n_page / N_page * 100, 0), " %)", sep='')) @@ -529,7 +546,10 @@ matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice # Fixes the height and width of the table according to # the number of station and the number of column to draw height = nsubCode - width = nbpMod * 2 * nPeriod_trend + nPeriod_trend + nPeriod_mean * nbpMod + nPeriod_mean + nbpMod + # width = nbpMod * 2 * nPeriod_trend + nPeriod_trend + nPeriod_mean * nbpMod + nPeriod_mean + nbpMod + + width = nbpMax * 2 * nPeriod_trend + nPeriod_trend + nPeriod_mean * nbpMax + nPeriod_mean + nbpMax + # Fixes the size of the plot area to keep proportion right options(repr.plot.width=width, repr.plot.height=height) @@ -544,7 +564,7 @@ matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice axis.ticks.y=element_blank(), axis.ticks.x=element_blank(), axis.title.y=element_blank(), - plot.margin=margin(t=5, r=5, b=5, l=5, unit="mm") + plot.margin=margin(t=0, r=0, b=0, l=0, unit="mm") ) # Extracts the name of the currently hydrological @@ -1093,8 +1113,59 @@ matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice dpi = 100 } + # If there is a foot note + if (foot_note) { + foot = foot_panel('tableau récapitulatif', + n_page, N_page, + resources_path, + AEAGlogo_file, INRAElogo_file, + FRlogo_file, foot_height) + + # Stores the map, the title and the colorbar in a list + P = list(mat, foot) + LM = matrix(c(1, + 2), + nrow=2, byrow=TRUE) + } else { + foot_height = 0 + # Stores the map, the title and the colorbar in a list + P = list(mat) + LM = matrix(c(1), + nrow=1, byrow=TRUE) + } + id_foot = 2 + + LMcol = ncol(LM) + LMrow = nrow(LM) + + LM = rbind(rep(99, times=LMcol), + LM, rep(99, times=LMcol)) + LMrow = nrow(LM) + LM = cbind(rep(99, times=LMrow), + LM, rep(99, times=LMrow)) + LMcol = ncol(LM) + + margin_height = 0.5 + + row_height = (height - 2*margin_height - foot_height) / (LMrow - 3) + + Hcut = LM[, 2] + heightLM = rep(row_height, times=LMrow) + heightLM[Hcut == id_foot] = foot_height + heightLM[Hcut == 99] = margin_height + + col_width = (width - 2*margin_height) / (LMcol - 2) + + Wcut = LM[(nrow(LM)-1),] + widthLM = rep(col_width, times=LMcol) + widthLM[Wcut == 99] = margin_height + + # Arranges the graphical object + plot = grid.arrange(grobs=P, layout_matrix=LM, + heights=heightLM, widths=widthLM) + # Saving - ggsave(plot=mat, + ggsave(plot=plot, path=outdirTmp, filename=paste(outnameTmp, '_', type, diff --git a/script.R b/script.R index e93ac8c4946b64c6f10a5be5c1d74d73435709f9..99f6c6777004e9f8e9a236a3fd8dcc506bd831e3 100644 --- a/script.R +++ b/script.R @@ -55,20 +55,19 @@ filedir = # Name of the file that will be analysed from the BH directory # (if 'all', all the file of the directory will be chosen) filename = - # "" - - c( + "" + # c( # "S2235610_HYDRO_QJM.txt", # "P1712910_HYDRO_QJM.txt", # "P0885010_HYDRO_QJM.txt", # "O5055010_HYDRO_QJM.txt", # "O0384010_HYDRO_QJM.txt", # "S4214010_HYDRO_QJM.txt", - "Q7002910_HYDRO_QJM.txt" + # "Q7002910_HYDRO_QJM.txt" # "O3035210_HYDRO_QJM.txt" # "O0554010_HYDRO_QJM.txt", # "O1584610_HYDRO_QJM.txt" - ) + # ) ## AGENCE EAU ADOUR GARONNE SELECTION @@ -78,8 +77,8 @@ AGlistdir = "" AGlistname = - "" - # "Liste-station_RRSE.docx" + # "" + "Liste-station_RRSE.docx" ## NIVALE SELECTION @@ -103,8 +102,8 @@ period1 = c("1968-01-01", "1988-12-31") period2 = c("2000-01-01", "2020-12-31") mean_period = list(period1, period2) -# p value thresold -p_thresold = 0.1 +# alpha the risk +alpha = 0.1 # Sampling span of the data sampleSpan = c('05-01', '11-30') @@ -192,8 +191,7 @@ if (AGlistname != ""){ 'longueur_serie')) # Get filenames of the selection - filename = df_selec_AG[df_selec_AG$ok,]$filename - + filename = df_selec_AG[df_selec_AG$ok,]$filename # Extract metadata about selected stations df_meta_AG = extract_meta(computer_data_path, filedir, filename) # Extract data about selected stations @@ -210,7 +208,6 @@ if (INlistname != ""){ # Get filenames of the selection filename = df_selec_IN[df_selec_IN$ok,]$filename - # Extract metadata about selected stations df_meta_IN = extract_meta(computer_data_path, filedir, filename) # Extract data about selected stations @@ -239,37 +236,37 @@ df_meta = get_lacune(df_data, df_meta) df_meta = get_hydrograph(df_data, df_meta, period=mean_period[[1]])$meta ### 3.2. Trend analysis -# QA trend -res_QAtrend = get_QAtrend(df_data, df_meta, - period=trend_period, - p_thresold=p_thresold) - -# QMNA tend -res_QMNAtrend = get_QMNAtrend(df_data, df_meta, - period=trend_period, - p_thresold=p_thresold, - sampleSpan=sampleSpan) - -# VCN10 trend -res_VCN10trend = get_VCN10trend(df_data, df_meta, - period=trend_period, - p_thresold=p_thresold, - sampleSpan=sampleSpan) - -# Start date for low water trend -res_DEBtrend = get_DEBtrend(df_data, df_meta, - period=trend_period, - p_thresold=p_thresold, - sampleSpan=sampleSpan, - thresold_type='VCN10', - select_longest=TRUE) -# res_DEBtrend = read_listofdf(resdir, 'res_DEBtrend') - -# Center date for low water trend -res_CENtrend = get_CENtrend(df_data, df_meta, - period=trend_period, - p_thresold=p_thresold, - sampleSpan=sampleSpan) +# # QA trend +# res_QAtrend = get_QAtrend(df_data, df_meta, +# period=trend_period, +# alpha=alpha) + +# # QMNA tend +# res_QMNAtrend = get_QMNAtrend(df_data, df_meta, +# period=trend_period, +# alpha=alpha, +# sampleSpan=sampleSpan) + +# # VCN10 trend +# res_VCN10trend = get_VCN10trend(df_data, df_meta, +# period=trend_period, +# alpha=alpha, +# sampleSpan=sampleSpan) + +# # Start date for low water trend +# res_DEBtrend = get_DEBtrend(df_data, df_meta, +# period=trend_period, +# alpha=alpha, +# sampleSpan=sampleSpan, +# thresold_type='VCN10', +# select_longest=TRUE) +# # res_DEBtrend = read_listofdf(resdir, 'res_DEBtrend') + +# # Center date for low water trend +# res_CENtrend = get_CENtrend(df_data, df_meta, +# period=trend_period, +# alpha=alpha, +# sampleSpan=sampleSpan) ### 3.3. Break analysis # df_break = get_break(res_QAtrend$data, df_meta) @@ -289,7 +286,7 @@ df_shapefile = ini_shapefile(computer_data_path, fr_shpdir, fr_shpname, bs_shpdir, bs_shpname, sbs_shpdir, sbs_shpname, - rv_shpdir, rv_shpname, riv=FALSE) + rv_shpdir, rv_shpname, riv=TRUE) ### 4.1. Simple time panel to criticize station data # Plot time panel of debit by stations @@ -307,7 +304,7 @@ df_shapefile = ini_shapefile(computer_data_path, ### 4.2. Analysis layout datasheet_layout(toplot=c( 'datasheet', - # 'matrix', + 'matrix', 'map' ), df_meta=df_meta,