matrix.R

# \\\
# Copyright 2021-2022 Louis Héraut*1
#
# *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/matrix.R
#
#


## 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) {

    # Number of variable/plot
    nbp = length(list_df2plot)
    
    # Get all different stations code
    Code = levels(factor(df_meta$code))
    nCode = length(Code)

    # Gets a trend example
    df_trend = list_df2plot[[1]]$trend

    # Convert 'trend_period' to list
    trend_period = as.list(trend_period)
    # Number of trend period
    nPeriod_trend = length(trend_period)

    # Fix the maximal number of period to the minimal possible
    nPeriod_max = 0
    # For all code 
    for (code in Code) {
        # 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))

        # If the number of period for the trend is greater
        # than the current max period, stocks it
        if (nPeriod > nPeriod_max) {
            nPeriod_max = nPeriod
        }
    }

    # 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 trend period
        for (i in 1:nPeriod_trend) {
            # Stocks period
            Periods = append(Periods, 
                             paste(Start[i],
                                   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
    for (j in 1:nPeriod_trend) {
        # 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
                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
                }
            }
        }
    }
    # Computes 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)

    # Blank vectors to store info about trend analyses
    Periods_trend = c()
    NPeriod_trend = c()
    Type_trend = list()
    Code_trend = c()
    Pthresold_trend = c()
    TrendMean_trend = c()
    DataMean_trend = c()
    Fill_trend = c()
    Color_trend = c()
    # For all the trend period
    for (j in 1:nPeriod_trend) {
        # For all code
        for (code in Code) {
            # 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
                # Extract the variable of the plot
                type = list_df2plot[[i]]$type
                # 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
                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){
                    # Gets the color associated to the mean trend
                    color_res = get_color(trendMean, 
                                          minTrendMean[j, i],
                                          maxTrendMean[j, i],
                                          palette_name='perso',
                                          reverse=TRUE)
                    # Specifies the color fill and contour of
                    # table cells
                    fill = color_res
                    color = 'white'
                    Pthresold = p_thresold
                # Otherwise it is not significative
                } else { 
                    fill = 'white'
                    color = 'grey85'  
                    Pthresold = NA
                }

                # Stores info needed to plot
                Periods_trend = append(Periods_trend, Periods)
                NPeriod_trend = append(NPeriod_trend, j)
                Type_trend = append(Type_trend, type)
                Code_trend = append(Code_trend, code)
                Pthresold_trend = append(Pthresold_trend, Pthresold)
                TrendMean_trend = append(TrendMean_trend, trendMean)
                DataMean_trend = append(DataMean_trend, dataMean)
                Fill_trend = append(Fill_trend, fill)
                Color_trend = append(Color_trend, color)
            }
        }
    }

    # If there is a 'mean_period'
    if (!is.null(mean_period)) {
        # Blank vectors to store info about mean analyses
        Periods_mean = c()
        NPeriod_mean = c()
        Type_mean = list()
        Code_mean = c()
        DataMean_mean = c()
        BreakMean_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
        BreakMean_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 mean period
        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
                    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)
                    # Creates a period name
                    Periods = paste(Datemin, Datemax,
                                    sep=' / ')

                    # Mean of the flow over the sub period
                    dataMean = mean(df_data_code_per$Qm3s,
                                    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
                    }
                    # Normalises the break by the mean of the
                    # initial period
                    BreakMean = Break / dataMeantmp[i, k]
                    # Stores the result
                    BreakMean_code[j, i, k] = BreakMean
                    # Stores temporarily the mean of the current period
                    dataMeantmp[i, k] = dataMean
                    
                    # Stores info needed to plot
                    Periods_mean = append(Periods_mean, Periods)
                    NPeriod_mean = append(NPeriod_mean, j)
                    Type_mean = append(Type_mean, type)
                    Code_mean = append(Code_mean, code)
                    DataMean_mean = append(DataMean_mean, dataMean)
                    BreakMean_mean = append(BreakMean_mean,
                                            BreakMean)
                }
            }
        }
        # Computes the min and the max of the mean trend for
        # all the station
        minBreakMean = apply(BreakMean_code, c(1, 2),
                             min, na.rm=TRUE)
        maxBreakMean = apply(BreakMean_code, c(1, 2),
                             max, na.rm=TRUE)
        # Blanks vector to store color info
        Fill_mean = c()
        Color_mean = c()
        # Index to count over all break computed
        ii = 1
        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 break mean
                    BreakMean = BreakMean_mean[ii]
                    # Gets the color associated
                    color_res = get_color(BreakMean, 
                                          minBreakMean[j, i],
                                          maxBreakMean[j, i],
                                          palette_name='perso',
                                          reverse=TRUE)
                    # Gets the fill and contour color
                    fill = color_res
                    color = 'white'
                    # Stores it
                    Fill_mean = append(Fill_mean, fill)
                    Color_mean = append(Color_mean, color)
                    # Passes to the next index
                    ii = ii + 1
                }
            }
        }
    }

    # If the slice option is not specified, the info for all
    # stations will be draw on the same page 
    if (is.null(slice)) {
        slice = nCode
    }

    # Extracts each possibilities of first letter of station code
    firstLetter = levels(factor(substr(Code, 1, 1)))
    # For all the available first letter
    for (fL in firstLetter) {
        # Print the matrix name
        print(paste('Matrix for region :', fL))

        # Get only station code with the same first letter 
        subCodefL = Code[substr(Code, 1, 1) == fL]
        # Counts the number of station in it
        nsubCodefL = length(subCodefL)
        # Computes the number of pages needed to plot all stations
        nMat = as.integer(nsubCodefL/slice) + 1
        # For all the pages
        for (imat in 1:nMat) {
            # Extracts the station for the current page
            subCode = subCodefL[(slice*(imat-1)+1):(slice*imat)]
            # Removes NA stations
            subCode = subCode[!is.na(subCode)]
            # Gets the number of station for the page
            nsubCode = length(subCode)

            # Creates logical vector to select only info about
            # stations that will be plot on the page
            CodefL_trend = Code_trend %in% subCode
            # Extracts those info
            subPeriods_trend = Periods_trend[CodefL_trend]
            subNPeriod_trend = NPeriod_trend[CodefL_trend]
            subType_trend = Type_trend[CodefL_trend]
            subCode_trend = Code_trend[CodefL_trend]
            subPthresold_trend = Pthresold_trend[CodefL_trend]
            subTrendMean_trend = TrendMean_trend[CodefL_trend]
            subDataMean_trend = DataMean_trend[CodefL_trend]
            subFill_trend = Fill_trend[CodefL_trend]
            subColor_trend = Color_trend[CodefL_trend]

            # Same for mean difference analysis
            CodefL_mean = Code_mean %in% subCode
            # Extracts right info
            subPeriods_mean = Periods_mean[CodefL_mean]
            subNPeriod_mean = NPeriod_mean[CodefL_mean]
            subType_mean = Type_mean[CodefL_mean]
            subCode_mean = Code_mean[CodefL_mean]
            subDataMean_mean = DataMean_mean[CodefL_mean]
            subBreakMean_mean = BreakMean_mean[CodefL_mean]
            subFill_mean = Fill_mean[CodefL_mean]
            subColor_mean = Color_mean[CodefL_mean]

            # Extracts the name of the currently hydrological
            # region plotted
            title = df_meta[df_meta$code == subCode[1],]$region_hydro
            
            ### Plot ###
            # Fixes the height and width of the table according to
            # the number of station and the number of column to draw
            height = nsubCode
            width = nbp * 2 * nPeriod_trend + nPeriod_trend + nPeriod_mean * nbp + nPeriod_mean + nbp        

            # Fixes the size of the plot area to keep proportion right
            options(repr.plot.width=width, repr.plot.height=height)

            # Open a new plot with a personalise theme
            mat = ggplot() + theme_ash + 
                # Modification of theme in order to remove axis
                theme(
                    panel.border=element_blank(),
                    axis.text.x=element_blank(),
                    axis.text.y=element_blank(),
                    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")
                )
            # Postion and name of the title
            xt = 1 - 6
            yt = height + 2
            Title = bquote(bold(.(title)))
            # Writes the title
            mat = mat +
                annotate("text", x=xt, y=yt,
                         label=Title,
                         hjust=0, vjust=1, 
                         size=6, color="#00A3A8")

            ### Trend ###
            # For all the trend period
            for (j in 1:nPeriod_trend) {
                # Extracts the info to plot associated to the
                # right period
                Type_trend_per =
                    subType_trend[subNPeriod_trend == j]
                Code_trend_per =
                    subCode_trend[subNPeriod_trend == j]
                Pthresold_trend_per =
                    subPthresold_trend[subNPeriod_trend == j]
                TrendMean_trend_per =
                    subTrendMean_trend[subNPeriod_trend == j]
                DataMean_trend_per =
                    subDataMean_trend[subNPeriod_trend == j]
                Fill_trend_per =
                    subFill_trend[subNPeriod_trend == j]
                Color_trend_per =
                    subColor_trend[subNPeriod_trend == j]

                # Converts the vector of hydrological variable to
                # a vector of integer associated to those variable
                Xtmp = as.integer(factor(as.character(Type_trend_per)))
                # Computes X position of the column of the period dates
                Xc = j + (j - 1)*nbp*2
                # Computes X positions of columns of the mean of
                # variables
                Xm = Xtmp + (j - 1)*nbp*2 + j
                # Computes X positions of columns of the mean trend
                X = Xtmp + (j - 1)*nbp*2 + nbp + j
                
                # Computes Y positions of each line for each station
                Y = as.integer(factor(Code_trend_per))

                # Position of a line to delimite periods
                x = Xc - 0.4
                xend = X[length(X)] + 0.25
                y = height + 1
                yend = height + 1
                # Drawing of the line
                mat = mat +
                    annotate("segment",
                             x=x, xend=xend,
                             y=y, yend=yend, 
                             color="grey40", size=0.35)

                # Position of the name of the current period
                yt = y + 0.15
                Start = trend_period[[j]][1]
                End = trend_period[[j]][2]
                # Name of the period
                periodName = bquote(bold('Période')~bold(.(as.character(j))))
                # Naming the period
                mat = mat +
                    annotate("text", x=x, y=yt,
                             label=periodName,
                             hjust=0, vjust=0.5, 
                             size=3, color='grey40')

                # For all the variable
                for (i in 1:length(X)) {
                    mat = mat +
                        # Plots circles for mean trends
                        gg_circle(r=0.45, xc=X[i], yc=Y[i],
                                  fill=Fill_trend_per[i],
                                  color=Color_trend_per[i]) +
                        # Plots circles for mean of variables
                        gg_circle(r=0.45, xc=Xm[i], yc=Y[i],
                                  fill='white', color='grey40') +
                        # Plots circles for the column of period dates
                        gg_circle(r=0.45, xc=Xc, yc=Y[i],
                                  fill='white', color='grey40') 
                }
                # For all the mean trend on this period
                for (i in 1:length(TrendMean_trend_per)) {
                    # Extracts the value of the mean trend
                    trendMean = TrendMean_trend_per[i]
                    # Converts it to the right format with two
                    # significant figures
                    trendMeanC = signif(trendMean*100, 2)

                    # If it is significative
                    if (!is.na(Pthresold_trend_per[i])) {
                        # The text color is white
                        Tcolor = 'white'
                    # Otherwise
                    } else {
                        # The text is grey
                        Tcolor = 'grey85'
                    }
                    
                    # Same for the mean of a variable over
                    # the current period
                    dataMean = DataMean_trend_per[i]
                    dataMeanC = signif(dataMean, 2)

                    mat = mat +
                        # Writes the mean trend
                        annotate('text', x=X[i], y=Y[i],
                                 label=trendMeanC,
                                 hjust=0.5, vjust=0.5, 
                                 size=3, color=Tcolor) + 
                        # Writes the mean of the associated variable
                        annotate('text', x=Xm[i], y=Y[i],
                                 label=dataMeanC,
                                 hjust=0.5, vjust=0.5, 
                                 size=3, color='grey40')
                }
                
                mat = mat +
                    # Writes a name for the period dates column
                    annotate('text', x=Xc, y=max(Y) + 0.85,
                             label=bquote(bold('Début')),
                             hjust=0.5, vjust=0.5, 
                             size=3, color='grey20') + 
                    annotate('text', x=Xc, y=max(Y) + 0.6,
                             label=bquote(bold('Fin')),
                             hjust=0.5, vjust=0.5, 
                             size=3, color='grey20')
                
                # For all variable
                for (i in 1:nbp) {
                    # Extract the variable of the plot
                    type = list_df2plot[[i]]$type
                    mat = mat +
                        # Writes the type of the variable
                        annotate('text', x=X[i], y=max(Y) + 0.82,
                                 label=bquote(.(type)),
                                 hjust=0.5, vjust=0.5, 
                                 size=3.25, color='grey20') +
                        # Writes the unit of the variable
                        annotate('text', x=X[i], y=max(Y) + 0.6,
                                 label=bquote('[%.'*ans^{-1}*']'),
                                 hjust=0.5, vjust=0.5, 
                                 size=2, color='grey40') +
                        # Writes the type of the mean variable
                        annotate('text', x=Xm[i], y=max(Y) + 0.82,
                                 label=bquote('µ'*.(type)),
                                 hjust=0.5, vjust=0.5, 
                                 size=3.25, color='grey20') +
                        # Writes the unit of the mean variable
                        annotate('text', x=Xm[i], y=max(Y) + 0.6,
                                 label=bquote('['*m^3*'.'*s^{-1}*']'),
                                 hjust=0.5, vjust=0.5, 
                                 size=2, color='grey40')
                }

                # For all the station on the page
                for (k in 1:nsubCode) {
                    # Gets the code
                    code = subCode[k]

                    label = Periods_trend[subNPeriod_trend == j
                                          & subCode_trend == code][1]
                    periodStart = substr(label, 1, 4)
                    periodEnd = substr(label, 14, 17)
                    
                    mat = mat +
                        annotate('text', x=Xc, y=k + 0.13,
                                 label=bquote(bold(.(periodStart))),
                                 hjust=0.5, vjust=0.5, 
                                 size=3, color='grey40') + 
                        
                        annotate('text', x=Xc, y=k - 0.13,
                                 label=bquote(bold(.(periodEnd))),
                                 hjust=0.5, vjust=0.5, 
                                 size=3, color='grey40')
                }
            }

            ### Mean ###
            for (j in 1:nPeriod_mean) {

                Type_mean_per =
                    subType_mean[subNPeriod_mean == j]
                Code_mean_per =
                    subCode_mean[subNPeriod_mean == j]
                DataMean_mean_per =
                    subDataMean_mean[subNPeriod_mean == j]
                BreakMean_mean_per =
                    subBreakMean_mean[subNPeriod_mean == j]
                Fill_mean_per =
                    subFill_mean[subNPeriod_mean == j]
                Color_mean_per =
                    subColor_mean[subNPeriod_mean == j]
                
                Xtmp_mean = as.integer(factor(as.character(Type_mean_per)))

                Xc_mean = j + (j - 1)*nbp + X[length(X)]
                Xm_mean = Xtmp_mean + (j - 1)*nbp + j + X[length(X)]
                Xr_mean = Xtmp_mean + (j - 1)*nbp*2 + j + X[length(X)]
                
                Y_mean = as.integer(factor(Code_mean_per))

                x = Xc_mean - 0.4
                xend = Xm_mean[length(Xm_mean)] + 0.25
                y = height + 1
                yend = height + 1

                mat = mat +
                    annotate("segment",
                             x=x, xend=xend,
                             y=y, yend=yend, 
                             color="grey40", size=0.35)
                
                yt = y + 0.15
                Start = mean_period[[j]][1]
                End = mean_period[[j]][2]
                
                periodName = bquote(bold('Période')~bold(.(as.character(j+nPeriod_trend))))

                mat = mat +
                    annotate("text", x=x, y=yt,
                             label=periodName,
                             hjust=0, vjust=0.5, 
                             size=3, color='grey40')
                
                if (j > 1) {
                    x = Xr_mean[1] - 0.4
                    xend = Xr_mean[length(Xr_mean)] + 0.25

                    mat = mat +
                        annotate("segment",
                                 x=x, xend=xend,
                                 y=y, yend=yend, 
                                 color="grey40", size=0.35)
                    
                    breakName =  bquote(bold('Écart')~bold(.(as.character(j-1+nPeriod_trend)))*bold('-')*bold(.(as.character(j+nPeriod_trend))))
                    
                    mat = mat +
                        annotate("text", x=x, y=yt,
                                 label=breakName,
                                 hjust=0, vjust=0.5, 
                                 size=3, color='grey40')
                }
                
                for (i in 1:length(Xm_mean)) {
                    mat = mat +
                        
                        gg_circle(r=0.45, xc=Xm_mean[i], yc=Y[i],
                                  fill='white', color='grey40') +
                        
                        gg_circle(r=0.45, xc=Xc_mean, yc=Y[i],
                                  fill='white', color='grey40')
                    
                    if (j > 1) {
                        mat = mat +
                            gg_circle(r=0.45, xc=Xr_mean[i], yc=Y[i],
                                      fill=Fill_mean_per[i],
                                      color=Color_mean_per[i])
                    }
                }

                for (i in 1:length(DataMean_mean_per)) {
                    dataMean = signif(DataMean_mean_per[i], 2)
                    
                    mat = mat +
                        annotate('text', x=Xm_mean[i], y=Y[i],
                                 label=dataMean,
                                 hjust=0.5, vjust=0.5, 
                                 size=3, color='grey40')

                    if (j > 1) {
                        BreakMean = BreakMean_mean_per[i]
                        BreakC = signif(BreakMean*100, 2)
                        
                        mat = mat +
                            annotate('text', x=Xr_mean[i], y=Y[i],
                                     label=BreakC,
                                     hjust=0.5, vjust=0.5, 
                                     size=3, color='white')   
                    }
                }
                
                mat = mat +
                    annotate('text', x=Xc_mean, y=max(Y) + 0.85,
                             label=bquote(bold('Début')),
                             hjust=0.5, vjust=0.5, 
                             size=3, color='grey20') + 
                    
                    annotate('text', x=Xc_mean, y=max(Y) + 0.6,
                             label=bquote(bold('Fin')),
                             hjust=0.5, vjust=0.5, 
                             size=3, color='grey20')
                # For all variable
                for (i in 1:nbp) {
                    # Extract the variable of the plot
                    type = list_df2plot[[i]]$type
                    mat = mat +
                        annotate('text', x=Xm_mean[i], y=max(Y) + 0.82,
                                 label=bquote('µ'*.(type)),
                                 hjust=0.5, vjust=0.5, 
                                 size=3.25, color='grey20') +
                        
                        annotate('text', x=Xm_mean[i], y=max(Y) + 0.6,
                                 label=bquote('['*m^3*'.'*s^{-1}*']'),
                                 hjust=0.5, vjust=0.5, 
                                 size=2, color='grey40')

                    if (j > 1) {
                        mat = mat +
                            annotate('text', x=Xr_mean[i], y=max(Y) + 0.82,
                                     label=bquote('d'*.(type)),
                                     hjust=0.5, vjust=0.5, 
                                     size=3.25, color='grey20') +
                            
                            annotate('text', x=Xr_mean[i], y=max(Y) + 0.6,
                                     label=bquote('[%]'),
                                     hjust=0.5, vjust=0.5, 
                                     size=2, color='grey40')
                            
                    }
                }

                for (k in 1:nsubCode) {
                    code = subCode[k]
                    label = Periods_mean[subNPeriod_mean == j
                                         & subCode_mean == code][1]
                    periodStart = substr(label, 1, 4)
                    periodEnd = substr(label, 14, 17)

                    mat = mat +
                        annotate('text', x=Xc_mean, y=k + 0.13,
                                 label=bquote(bold(.(periodStart))),
                                 hjust=0.5, vjust=0.5, 
                                 size=3, color='grey40') + 
                        
                        annotate('text', x=Xc_mean, y=k - 0.13,
                                 label=bquote(bold(.(periodEnd))),
                                 hjust=0.5, vjust=0.5, 
                                 size=3, color='grey40')
                }            
            }
            
            ### Code ###
            for (k in 1:nsubCode) {
                # Gets the code
                code = subCode[k]
                name = df_meta[df_meta$code == code,]$nom
                ncharMax = 38
                if (nchar(name) > ncharMax) {
                    name = paste(substr(name, 1, ncharMax), '...', sep='')
                }

                mat = mat +
                    annotate('text', x=0.3, y=k + 0.14,
                             label=bquote(bold(.(code))),
                             hjust=1, vjust=0.5, 
                             size=3.5, color="#00A3A8") +
                    
                    annotate('text', x=0.3, y=k - 0.14,
                             label=name,
                             hjust=1, vjust=0.5, 
                             size=3.5, color="#00A3A8")
            }

            ### Environment ###
            mat = mat +
                # Fixed coordinate system
                coord_fixed() +
                
                scale_x_continuous(limits=c(1 - rel(6), 
                                            width + rel(0.5)),
                                   expand=c(0, 0)) + 
                
                scale_y_continuous(limits=c(1 - rel(0.5), 
                                            height + rel(2)),
                                   expand=c(0, 0))
            
            # Saving matrix plot
            if (A3) {
                width = 42
                height = 29.7
                dpi = 300
            } else {
                width = 29.7
                height = 21
                dpi = 100
            }
            
            ggsave(plot=mat, 
                   path=outdirTmp,
                   filename=paste(outnameTmp, '_', fL, imat, '.pdf',
                                  sep=''),
                   width=width, height=height, units='cm', dpi=dpi)            
        }
    }           
}