diff --git a/plotting/map.R b/plotting/map.R index a7ccdbd36169662bf3e11c1814433c172d33ccd3..c918ea0e5123385faf2574399077a36fae28edc5 100644 --- a/plotting/map.R +++ b/plotting/map.R @@ -27,6 +27,7 @@ ## 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, verbose=TRUE) { # Extract shapefiles @@ -100,14 +101,15 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='' substr(End[i], 1, 4), sep=' / ')) } - + # Stores time info by station Start_code[[j]] = Start End_code[[j]] = End Code_code[[j]] = code Periods_code[[j]] = Periods - } - + + # Blank array to store mean of the trend for each + # station, perdiod and variable TrendMean_code = array(rep(1, nPeriod_max*nbp*nCode), dim=c(nPeriod_max, nbp, nCode)) # For all the period @@ -130,54 +132,72 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='' # Extracts the trend corresponding to the code df_trend_code = df_trend[df_trend$code == code,] + # Gets the associated time info Start = Start_code[Code_code == code][[1]][j] End = End_code[Code_code == code][[1]][j] Periods = Periods_code[Code_code == code][[1]][j] + # 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,] } - + + # Computes the mean of the data on the period dataMean = mean(df_data_code_per$Qm3s, na.rm=TRUE) + # Normalises the trend value by the mean of the data trendMean = df_trend_code_per$trend / dataMean + # If the p value is under the threshold if (df_trend_code_per$p <= p_threshold){ + # Stores the mean trend TrendMean_code[j, i, k] = trendMean + # Otherwise } else { + # Do not stocks it TrendMean_code[j, i, k] = NA } } } } - + # Compute the min and the max of the mean trend for all the station minTrendMean = apply(TrendMean_code, c(1, 2), min, na.rm=TRUE) maxTrendMean = apply(TrendMean_code, c(1, 2), max, na.rm=TRUE) - ncolor = 256 + # 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 } - # Extract the variable of the plot + # Extracts the variable of the plot type = list_df2plot[[i]]$type + # Createsa name for the map outname = paste('map_', type, sep='') + # If there is the verbose option if (verbose) { + # Prints the name of the map print(paste('Map for variable :', type)) } + # 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 @@ -186,19 +206,20 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='' sizebs = 0.3 sizerv = 0.2 } - + + # Open a new plot with the personalise theme map = ggplot() + theme_void() + - # theme(plot.background=element_rect(fill=NA, # color="#EC4899")) + # Fixed coordinate system coord_fixed() + - + # 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), @@ -206,30 +227,44 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='' } 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 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(280000, 790000) ylim = c(6110000, 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(7, 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(5, xlim, shift=TRUE) xint = c(0, 100*1E3) ymin = gpct(1, ylim, shift=TRUE) @@ -239,38 +274,47 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='' } 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) { 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 + + 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)) { - map = map + + # Sets all margins to 0 + map = map + theme(plot.margin=margin(t=0, r=0, b=0, l=0, unit="mm")) + # Otherwise } else { - map = map + + # 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() @@ -314,13 +358,13 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='' maxTrendMean[idPer, i], palette_name='perso', reverse=TRUE, - ncolor=ncolor) + ncolor=256) palette_res = get_palette(minTrendMean[idPer, i], maxTrendMean[idPer, i], palette_name='perso', reverse=TRUE, - ncolor=ncolor, + ncolor=256, nbTick=nbTick) if (df_trend_code_per$p <= p_threshold){ @@ -355,12 +399,14 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='' 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 + 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 { plot_map_codeNo = plot_map[plot_map$code != codeLight,] plot_map_code = plot_map[plot_map$code == codeLight,] @@ -513,7 +559,6 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='' yTrend = c(yTrend, rep(mids[ii], times=counts[ii])) } - ## No touch distribution ## # start_hist = 1.25 # min_xsep = 0.15 @@ -533,10 +578,11 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='' # } # } - + # 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, @@ -545,11 +591,12 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='' pal = pal + + # Arrow to show a worsening of the situation geom_segment(aes(x=2.7, y=valNorm*0.75, xend=2.7, yend=valNorm*0.25), color='grey50', size=0.3, arrow=arrow(length=unit(2, "mm"))) + - + # Text associated to the arrow annotate('text', x=2.8, y=valNorm*0.5, label= "Plus sévère", @@ -559,23 +606,24 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='' pal = pal + - + # X axis of the colorbar scale_x_continuous(limits=c(-1, 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=5, b=5, l=0, unit="mm")) - + # Stores the map, the title and the colorbarin a list Map = list(map, title, pal) - + # Arranges the graphical object plot = grid.arrange(grobs=Map, layout_matrix= matrix(c(1, 1, 1, 2, 1, 1, 1, 3), nrow=2, byrow=TRUE)) + # If there is no specified station code to highlight (mini map) if (is.null(codeLight)) { # Saving matrix plot ggsave(plot=plot, @@ -584,6 +632,7 @@ map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='' width=29.7, height=21, units='cm', dpi=100) } } + # Returns the map object return (map) } diff --git a/plotting/matrix.R b/plotting/matrix.R index 16750b96509ab76efcca743382ffdc29a0a5ff7c..ba442bb566d9df2558a1fea7b4d50ade401d8c94 100644 --- a/plotting/matrix.R +++ b/plotting/matrix.R @@ -101,14 +101,15 @@ matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice End[i], sep=' / ')) } - + # Stores time info by station Start_code[[j]] = Start End_code[[j]] = End Code_code[[j]] = code Periods_code[[j]] = Periods - } - + + # Blank array to store mean of the trend for each + # station, perdiod and variable TrendMean_code = array(rep(1, nPeriod_trend*nbp*nCode), dim=c(nPeriod_trend, nbp, nCode)) # For all the trend period @@ -131,35 +132,46 @@ matrix_panel = function (list_df2plot, df_meta, trend_period, mean_period, slice # Extracts the trend corresponding to the code df_trend_code = df_trend[df_trend$code == code,] + # Gets the associated time info Start = Start_code[Code_code == code][[1]][j] End = End_code[Code_code == code][[1]][j] Periods = Periods_code[Code_code == code][[1]][j] + # 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,] } - + + # Computes the mean of the data on the period dataMean = mean(df_data_code_per$Qm3s, na.rm=TRUE) + # Normalises the trend value by the mean of the data trendMean = df_trend_code_per$trend / dataMean + # If the p value is under the threshold if (df_trend_code_per$p <= p_threshold){ + # Stores the mean trend TrendMean_code[j, i, k] = trendMean + # Otherwise } else { + # Do not stocks it TrendMean_code[j, i, k] = NA } } } } - + # Compute the min and the max of the mean trend for all the station minTrendMean = apply(TrendMean_code, c(1, 2), min, na.rm=TRUE) maxTrendMean = apply(TrendMean_code, c(1, 2), max, na.rm=TRUE)