map.R 37 KB
Newer Older
Heraut Louis's avatar
Heraut Louis committed
# \\\
Heraut Louis's avatar
Heraut Louis committed
# Copyright 2021-2022 Louis Héraut*1
Heraut Louis's avatar
Heraut Louis committed
#
# *1   INRAE, France
#      louis.heraut@inrae.fr
#
# This file is part of ash R toolbox.
#
# ash R toolbox is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# ash R toolbox is distributed in the hope that it will be useful, but 
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with ash R toolbox.  If not, see <https://www.gnu.org/licenses/>.
# ///
#
#
# plotting/map.R
#
Heraut Louis's avatar
Heraut Louis committed
# Deals with the creation of a map for presenting the trend analysis of hydrological variables
Heraut Louis's avatar
Heraut Louis committed


Heraut Louis's avatar
Heraut Louis committed
## 1. MAP PANEL
Heraut Louis's avatar
Heraut Louis committed
# Generates a map plot of the tendancy of a hydrological variable
Heraut Louis's avatar
Heraut Louis committed
map_panel = function (list_df2plot, df_meta, df_shapefile, idPer_trend=1,
                      mean_period=NULL, outdirTmp='', codeLight=NULL,
                      margin=NULL, showSea=TRUE,
Heraut Louis's avatar
Heraut Louis committed
                      foot_note=FALSE,
                      foot_height=0, resources_path=NULL,
                      AEAGlogo_file=NULL, INRAElogo_file=NULL,
                      FRlogo_file=NULL,
                      verbose=TRUE) {
Heraut Louis's avatar
Heraut Louis committed

Heraut Louis's avatar
Heraut Louis committed
    # Extract shapefiles
Heraut Louis's avatar
Heraut Louis committed
    df_france = df_shapefile$france
    df_bassin = df_shapefile$bassin
Heraut Louis's avatar
Heraut Louis committed
    df_subbassin = df_shapefile$subbassin
Heraut Louis's avatar
Heraut Louis committed
    df_river = df_shapefile$river

Heraut Louis's avatar
Heraut Louis committed
    # Number of variable/plot
Heraut Louis's avatar
Heraut Louis committed
    nbp = length(list_df2plot)
  
    # Get all different stations code
    Code = levels(factor(df_meta$code))
    nCode = length(Code)
Heraut Louis's avatar
Heraut Louis committed

    # Gets a trend example
Heraut Louis's avatar
Heraut Louis committed
    df_trend = list_df2plot[[1]]$trend
    
    nPeriod_max = 0
    for (code in Code) {
Heraut Louis's avatar
Heraut Louis committed
        # Extracts the trend corresponding to the code
Heraut Louis's avatar
Heraut Louis committed
        df_trend_code = df_trend[df_trend$code == code,]
        
Heraut Louis's avatar
Heraut Louis committed
        # Extract start and end of trend periods
Heraut Louis's avatar
Heraut Louis committed
        Start = df_trend_code$period_start
        End = df_trend_code$period_end
Heraut Louis's avatar
Heraut Louis committed
        # Get the name of the different period
        UStart = levels(factor(Start))        
Heraut Louis's avatar
Heraut Louis committed
        UEnd = levels(factor(End))
Heraut Louis's avatar
Heraut Louis committed

        # Compute the max of different start and end
        # so the number of different period
Heraut Louis's avatar
Heraut Louis committed
        nPeriod = max(length(UStart), length(UEnd))
Heraut Louis's avatar
Heraut Louis committed

        # If the number of period for the trend is greater
        # than the current max period, stocks it
Heraut Louis's avatar
Heraut Louis committed
        if (nPeriod > nPeriod_max) {
            nPeriod_max = nPeriod
        }
    }
Heraut Louis's avatar
Heraut Louis committed
    # Blank array to store time info
Heraut Louis's avatar
Heraut Louis committed
    tab_Start =  array(rep('', nCode*nbp*nPeriod_max),
                       dim=c(nCode, nbp, nPeriod_max))
    tab_End = array(rep('', nCode*nbp*nPeriod_max),
                    dim=c(nCode, nbp, nPeriod_max))
    tab_Code = array(rep('', nCode*nbp*nPeriod_max),
                     dim=c(nCode, nbp, nPeriod_max))
    tab_Periods = array(rep('', nCode*nbp*nPeriod_max),
                        dim=c(nCode, nbp, nPeriod_max))
    
    # For all code
    for (k in 1:nCode) {
Heraut Louis's avatar
Heraut Louis committed
        # Gets the code
Heraut Louis's avatar
Heraut Louis committed
        code = Code[k]
Heraut Louis's avatar
Heraut Louis committed
        # For all the variable
Heraut Louis's avatar
Heraut Louis committed
        for (i in 1:nbp) {
            df_trend = list_df2plot[[i]]$trend
            # 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))
Heraut Louis's avatar
Heraut Louis committed

Heraut Louis's avatar
Heraut Louis committed
            # For all the period
            for (j in 1:nPeriod_max) {
                # Stocks period
                Periods = paste(Start[j],
                                End[j],
                                sep=' / ')
Heraut Louis's avatar
Heraut Louis committed

                # Saves the time info
Heraut Louis's avatar
Heraut Louis committed
                tab_Start[k, i, j] = as.character(Start[j])
                tab_End[k, i, j] = as.character(End[j])
                tab_Code[k, i, j] = code
                tab_Periods[k, i, j] = Periods
                
            }
Heraut Louis's avatar
Heraut Louis committed
        }
    }
Heraut Louis's avatar
Heraut Louis committed
    
Heraut Louis's avatar
Heraut Louis committed
    # Blank array to store mean of the trend for each
    # station, perdiod and variable
Heraut Louis's avatar
Heraut Louis committed
    TrendValue_code = array(rep(1, nPeriod_max*nbp*nCode),
Heraut Louis's avatar
Heraut Louis committed
                           dim=c(nPeriod_max, nbp, nCode))
Heraut Louis's avatar
Heraut Louis committed
    # For all the period
Heraut Louis's avatar
Heraut Louis committed
    for (j in 1:nPeriod_max) {
Heraut Louis's avatar
Heraut Louis committed
        # For all the code
Heraut Louis's avatar
Heraut Louis committed
        for (k in 1:nCode) {
Heraut Louis's avatar
Heraut Louis committed
            # Gets the code
Heraut Louis's avatar
Heraut Louis committed
            code = Code[k]
Heraut Louis's avatar
Heraut Louis committed
            # For all variable
Heraut Louis's avatar
Heraut Louis committed
            for (i in 1:nbp) {
Heraut Louis's avatar
Heraut Louis committed
                # Extracts the data corresponding to the
                # current variable
Heraut Louis's avatar
Heraut Louis committed
                df_data = list_df2plot[[i]]$data
Heraut Louis's avatar
Heraut Louis committed
                # Extracts the trend corresponding to the
                # current variable
Heraut Louis's avatar
Heraut Louis committed
                df_trend = list_df2plot[[i]]$trend
Heraut Louis's avatar
Heraut Louis committed
                # Extracts the type of the variable
                type = list_df2plot[[i]]$type
Heraut Louis's avatar
Heraut Louis committed
                alpha = list_df2plot[[i]]$alpha
Heraut Louis's avatar
Heraut Louis committed
                # Extracts the data corresponding to the code
                df_data_code = df_data[df_data$code == code,]
                # Extracts the trend corresponding to the code
Heraut Louis's avatar
Heraut Louis committed
                df_trend_code = df_trend[df_trend$code == code,]

Heraut Louis's avatar
Heraut Louis committed
                # Gets the associated time info
Heraut Louis's avatar
Heraut Louis committed
                Start = tab_Start[k, i, j]
                End = tab_End[k, i, j]
                Periods = tab_Periods[k, i, j]
Heraut Louis's avatar
Heraut Louis committed

Heraut Louis's avatar
Heraut Louis committed
                # Extracts the corresponding data for the period
Heraut Louis's avatar
Heraut Louis committed
                df_data_code_per =
                    df_data_code[df_data_code$Date >= Start 
                                 & df_data_code$Date <= End,]
Heraut Louis's avatar
Heraut Louis committed
                # Same for trend
Heraut Louis's avatar
Heraut Louis committed
                df_trend_code_per = 
                    df_trend_code[df_trend_code$period_start == Start 
                                  & df_trend_code$period_end == End,]

Heraut Louis's avatar
Heraut Louis committed
                # Computes the number of trend analysis selected
Heraut Louis's avatar
Heraut Louis committed
                Ntrend = nrow(df_trend_code_per)
Heraut Louis's avatar
Heraut Louis committed
                # If there is more than one trend on the same period
Heraut Louis's avatar
Heraut Louis committed
                if (Ntrend > 1) {
Heraut Louis's avatar
Heraut Louis committed
                    # Takes only the first because they are similar
Heraut Louis's avatar
Heraut Louis committed
                    df_trend_code_per = df_trend_code_per[1,]
                }
Heraut Louis's avatar
Heraut Louis committed
                # 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
                } else if (type == 'saisonnalité') {
                    trendValue = df_trend_code_per$trend
                }
Heraut Louis's avatar
Heraut Louis committed

Heraut Louis's avatar
Heraut Louis committed
                # If the p value is under the threshold
Heraut Louis's avatar
Heraut Louis committed
                if (df_trend_code_per$p <= alpha){
Heraut Louis's avatar
Heraut Louis committed
                    # Stores the mean trend
Heraut Louis's avatar
Heraut Louis committed
                    TrendValue_code[j, i, k] = trendValue
Heraut Louis's avatar
Heraut Louis committed
                # Otherwise
Heraut Louis's avatar
Heraut Louis committed
                } else {
Heraut Louis's avatar
Heraut Louis committed
                    # Do not stocks it
Heraut Louis's avatar
Heraut Louis committed
                    TrendValue_code[j, i, k] = NA
Heraut Louis's avatar
Heraut Louis committed
                }
            }
        }
    }
Heraut Louis's avatar
Heraut Louis committed
    # Compute the min and the max of the mean trend for all the station
Heraut Louis's avatar
Heraut Louis committed
    minTrendValue = apply(TrendValue_code, c(1, 2), min, na.rm=TRUE)
    maxTrendValue = apply(TrendValue_code, c(1, 2), max, na.rm=TRUE)    
Heraut Louis's avatar
Heraut Louis committed

Heraut Louis's avatar
Heraut Louis committed




    # 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
    }
    
Heraut Louis's avatar
Heraut Louis committed
    # Number of ticks for the colorbar
Heraut Louis's avatar
Heraut Louis committed
    nbTick = 10
Heraut Louis's avatar
Heraut Louis committed
    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)   
            }
            
Heraut Louis's avatar
Heraut Louis committed
            map = map +
Heraut Louis's avatar
Heraut Louis committed
                # 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)
            }
Heraut Louis's avatar
Heraut Louis committed

Heraut Louis's avatar
Heraut Louis committed
            # 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
            }
            
Heraut Louis's avatar
Heraut Louis committed
            map = map +
Heraut Louis's avatar
Heraut Louis committed
                # 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
Heraut Louis's avatar
Heraut Louis committed
                annotate("text",
Heraut Louis's avatar
Heraut Louis committed
                         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)
            }
Heraut Louis's avatar
Heraut Louis committed
            # 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 
                }
Heraut Louis's avatar
Heraut Louis committed
                # 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)
Heraut Louis's avatar
Heraut Louis committed
            }
Heraut Louis's avatar
Heraut Louis committed
            # 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
Heraut Louis's avatar
Heraut Louis committed
            # If it is a flow variable
            if (type == 'sévérité') {
Heraut Louis's avatar
Heraut Louis committed
                # Formatting of label in pourcent
                labTick = as.character(round(labTick*100, 2))
                # If it is a date variable
Heraut Louis's avatar
Heraut Louis committed
            } else if (type == 'saisonnalité') {
Heraut Louis's avatar
Heraut Louis committed
                # Formatting of label
                labTick = as.character(round(labTick, 2))
Heraut Louis's avatar
Heraut Louis committed
            
            # 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)"))
Heraut Louis's avatar
Heraut Louis committed
                }
Heraut Louis's avatar
Heraut Louis committed
            } else {
Heraut Louis's avatar
Heraut Louis committed
                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')
Heraut Louis's avatar
Heraut Louis committed
            }

Heraut Louis's avatar
Heraut Louis committed
            yUp = -20
            yNone = -29
Heraut Louis's avatar
Heraut Louis committed

Heraut Louis's avatar
Heraut Louis committed
            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
            }
            
Heraut Louis's avatar
Heraut Louis committed
            pal = pal +
Heraut Louis's avatar
Heraut Louis committed
                # 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,
Heraut Louis's avatar
Heraut Louis committed
                         size=3, color='grey40')

Heraut Louis's avatar
Heraut Louis committed
            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')
Heraut Louis's avatar
Heraut Louis committed
            }
Heraut Louis's avatar
Heraut Louis committed
            
            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)
Heraut Louis's avatar
Heraut Louis committed

Heraut Louis's avatar
Heraut Louis committed
            if (type == 'sévérité') {
                labelArrow = 'Plus sévère'
            } else if (type == 'saisonnalité') {
                labelArrow = 'Plus tôt'
            }
Heraut Louis's avatar
Heraut Louis committed

Heraut Louis's avatar
Heraut Louis committed
            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"))
Heraut Louis's avatar
Heraut Louis committed

Heraut Louis's avatar
Heraut Louis committed
            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
Heraut Louis's avatar
Heraut Louis committed

Heraut Louis's avatar
Heraut Louis committed
            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)
            }
Heraut Louis's avatar
Heraut Louis committed
        }
    }
Heraut Louis's avatar
Heraut Louis committed
    # Returns the map object
Heraut Louis's avatar
Heraut Louis committed
    return (map)
}