Skip to content
GitLab
Select Git revision
  • master default
Find file BlameHistoryPermalink
map.R 25.59 KiB
1 
# \\\
2 
# Copyright 2021-2022 Louis Héraut*1
3 
#
4 
# *1   INRAE, France
5 
#      louis.heraut@inrae.fr
6 
#
7 
# This file is part of ash R toolbox.
8 
#
9 
# ash R toolbox is free software: you can redistribute it and/or modify
10 
# it under the terms of the GNU General Public License as published by
11 
# the Free Software Foundation, either version 3 of the License, or (at
12 
# your option) any later version.
13 
#
14 
# ash R toolbox is distributed in the hope that it will be useful, but
15 
# WITHOUT ANY WARRANTY; without even the implied warranty of
16 
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17 
# General Public License for more details.
18 
#
19 
# You should have received a copy of the GNU General Public License
20 
# along with ash R toolbox.  If not, see <https://www.gnu.org/licenses/>.
21 
# ///
22 
#
23 
#
24 
# plotting/map.R
25 
#
26 
# Deals with the creation of a map for presenting the trend analysis of hydrological variables
27
28
29 
## 1. MAP PANEL
30 
# Generates a map plot of the tendancy of a hydrological variable
31 
map_panel = function (list_df2plot, df_meta, df_shapefile, idPer=1, outdirTmp='', codeLight=NULL, margin=NULL, showSea=TRUE, verbose=TRUE) {
32
33 
    # Extract shapefiles
34 
    df_france = df_shapefile$france
35 
    df_bassin = df_shapefile$bassin
36 
    df_river = df_shapefile$river
37
38 
    # Number of variable/plot
39 
    nbp = length(list_df2plot)
40
41 
    # Get all different stations code
42 
    Code = levels(factor(df_meta$code))
43 
    nCode = length(Code)
44
45 
    # Gets a trend example
46 
    df_trend = list_df2plot[[1]]$trend
47
48 
    nPeriod_max = 0
49 
    for (code in Code) {
50 
        # Extracts the trend corresponding to the code
51 
        df_trend_code = df_trend[df_trend$code == code,]
52
53 
        # Extract start and end of trend periods
54 
        Start = df_trend_code$period_start
55 
        End = df_trend_code$period_end
56 
        # Get the name of the different period
57 
        UStart = levels(factor(Start))
58 
        UEnd = levels(factor(End))
59
60 
        # Compute the max of different start and end
61 
        # so the number of different period
62 
        nPeriod = max(length(UStart), length(UEnd))
63
64 
        # If the number of period for the trend is greater
65 
        # than the current max period, stocks it
66 
        if (nPeriod > nPeriod_max) {
67 
            nPeriod_max = nPeriod
68 
        }
69 
    }
70
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
# Blank list to store time info by station code Start_code = vector(mode='list', length=nCode) End_code = vector(mode='list', length=nCode) Code_code = vector(mode='list', length=nCode) Periods_code = vector(mode='list', length=nCode) # For all the code for (j in 1:nCode) { # Gets the code code = Code[j] # Extracts the trend corresponding to the code df_trend_code = df_trend[df_trend$code == code,] # Extract start and end of trend periods Start = df_trend_code$period_start End = df_trend_code$period_end # Get the name of the different period UStart = levels(factor(Start)) UEnd = levels(factor(End)) # Compute the max of different start and end # so the number of different period nPeriod = max(length(UStart), length(UEnd)) # Vector to store trend period Periods = c() # For all the period for (i in 1:nPeriod_max) { # Stocks period Periods = append(Periods, paste(substr(Start[i], 1, 4), 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 for (j in 1:nPeriod_max) { # 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 # 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,] # 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
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
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) # 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 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 } else { sizefr = 0.35 sizebs = 0.3 sizerv = 0.2 } # Stores the coordonate system
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
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 + # 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) ymax = ymin + gpct(3, ylim) size = 2 sizekm = 1.5 }
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
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 + # 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 (code in Code) { # 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,] # Gets the associated time info Start = Start_code[Code_code == code][[1]][idPer] End = End_code[Code_code == code][[1]][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
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
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 # Computes the color associated to the mean trend color_res = get_color(trendMean, minTrendMean[idPer, i], maxTrendMean[idPer, i], palette_name='perso', reverse=TRUE, ncolor=256) # Computes the colorbar info palette_res = get_palette(minTrendMean[idPer, i], maxTrendMean[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 (trendMean >= 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 } # 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, trendMean) # 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
421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490
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='grey70', fill='grey70') + # 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='grey40', fill='grey40') } # 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 # Formatting of label in pourcent labTick = as.character(round(labTick*100, 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.3), y=c(0.05, 0.05)), size=0.6, color="#00A3A8") + # Writes title geom_shadowtext(data=tibble(x=-0.3, y=0.2, label=type), 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(-1, 1 + 3), expand=c(0, 0)) + # Y axis
491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560
scale_y_continuous(limits=c(0, 10), expand=c(0, 0)) + # Margin theme(plot.margin=margin(t=5, r=5, 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) pal = pal + # Name of the colorbar annotate('text', x=-0.3, y= valNorm + 23, label="Tendance", hjust=0, vjust=0.5, size=6, color='grey40') + # Unit legend of the colorbar annotate('text', x=-0.2, y= valNorm + 13, label=bquote(bold("% par an")), hjust=0, vjust=0.5, size=4, color='grey40') # For all the ticks for (j in 1:nbTick) { 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, 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%")),
561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630
hjust=0, vjust=0.5, size=3, color='grey40') # Normalises all the trend values for each station # according to the colorbar yTrend = (trend - minTrendMean[idPer, i]) / (maxTrendMean[idPer, i] - minTrendMean[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)) } # 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])) } ## No touch distribution ## # start_hist = 1.25 # min_xsep = 0.15 # min_ysep = 4 # xTrend = rep(start_hist, times=length(yTrend)) # for (ii in 1:length(yTrend)) { # yTrendtmp = yTrend # yTrendtmp[ii] = 1E99 # isinf_ysep = abs(yTrendtmp - yTrend[ii]) < min_ysep # if (any(isinf_ysep) & !all(xTrend[which(isinf_ysep)] > start_hist)) { # xTrend[ii] = max(xTrend[which(isinf_ysep)]) + min_xsep # } # } # Makes a tibble to plot the trend distribution plot_trend = tibble(xTrend=xTrend, yTrend=yTrend) pal = pal +
631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685
# 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) 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", angle=90, hjust=0.5, vjust=1, size=3, color='grey50') 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, path=outdirTmp, filename=paste(outname, '.pdf', sep=''), width=29.7, height=21, units='cm', dpi=100) } } # Returns the map object return (map) }